Morphology of Ion-Sputtered Surfaces 
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We derive a stochastic nonlinear continuum theory to describe the morphological evolution of 
amorphous surfaces eroded by ion bombardment. Starting from Sigmund's theory of sputter erosion, 
we calculate the coefficients appearing in the continuum equation in terms of the physical parame- 
ters characterizing the sputtering process. We analyze the morphological features predicted by the 
continuum theory, comparing them with the experimentally reported morphologies. We show that 
for short time scales, where the effect of nonlinear terms is negligible, the continuum theory predicts 
ripple formation. We demonstrate that in addition to relaxation by thermal surface diffusion, the 
sputtering process can also contribute to the smoothing mechanisms shaping the surface morphology. 
We explicitly calculate an effective surface diffusion constant characterizing this smoothing effect, 
and show that it is responsible for the low temperature ripple formation observed in various exper- 
iments. At long time scales the nonlinear terms dominate the evolution of the surface morphology. 
The nonlinear terms lead to the stabilization of the ripple wavelength and we show that, depending 
on the experimental parameters such as angle of incidence and ion energy, different morphologies 
can be observed: asymptotically, sputter eroded surfaces could undergo kinetic roughening, or can 
display novel ordered structures with rotated ripples. Finally, we discuss in detail the existing ex- 
perimental support for the proposed theory, and uncover novel features of the surface morphology 
and evolution, that could be directly tested experimentally. 
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PACS numbers; 79.20.Rf, 64.60.Ht, 68.35.Rh 



I. INTRODUCTION 

Sputtering is the removal of material from the surface 
of solids through the impact of energetic particles iQ-^. 
It is a widespread experimental technique, used in a large 
number of applications with a remarkable level of sophis- 
tication. It is a basic tool in surface analysis, depth profil- 
ing, sputter cleaning, micromachining, and sputter depo- 
sition. Perhaps the largest community of users is in the 
thin film and semiconductor fabrication areas, sputter 
erosion being routinely used for etching patterns impor- 
tant to the production of integrated circuits and device 
packaging. 

To have a better control over this important tool, we 
need to understand the effect of the sputtering process 
on the surface morphology. In many cases sputtering is 
routinely used to smooth out surface features. On the 
other hand, some investigations indicate that sputtering 
can also roughen the surface. Consequently, sputter ero- 
sion may have different effects on the surface, depending 
on many factors, such as incident ion energy, mass, angle 
of incidence, sputtered substrate temperature and mate- 
rial composition. The experimental results on the effect 
of sputter erosion on the surface morphology can be clas- 
sified in two main classes. There exists ample experimen- 
tal evidence that ion sputtering can lead to the develop- 
ment of periodic ripples on the surface P-pTj . Depending 
on the sputtered substrate and the sputtering conditions 
these ripples can be surprisingly straight and ordered. 
However, a number of recent investigations pl-pfl] have 



provided rather detailed and convincing experimental ev- 
idence that under certain experimental conditions ion 
eroded surfaces become rough, and the roughness follows 
the predictions of various scaling theories | |37|] . Moreover, 
these investigations did not find any evidence of ripple 
formation on the surface. Up to recently these two mor- 
phological features were treated separately and no unified 
theoretical framework describing these morphologies was 
available. 

The first widely accepted theoretical approach describ- 
ing the process of ripple formation on amorphous sub- 
strates was developed by Bradley and Harper (BH) Q . 
This theory is rather successful in predicting the ripple 
wavelength and orientation in agreement with numerous 
experimental observations. However, a number of exper- 
imental results have systematically eluded this theory. 
For example, the BH theory predicts an unlimited ex- 
ponential increase in ripple amplitude in contrast with 
the observed saturation of the surface width. Similarly, 
it cannot account for surface roughening, or for ripple 
orientations different from those defined by the incom- 
ing ion direction or perpendicular to it. Finally, recent 
experiments [^2|,^ have observed ripples whose wave- 
length is independent of the substrate temperature, and 
linear in the ion energy, in contrast with the BH predic- 
tion of a ripple wavelength which depends exponentially 
with temperature and decreases with ion energy. 

In the light of the accumulated experimental results, 
it is clear that a theory going beyond the BH approach 
is required, motivating the results described in this pa- 
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per. Thus here we investigate the morphology of ion- 
sputtered amorphous surfaces aiming to describe in an 
unified framework the dynamic and scahng behavior of 
the experimentally observed surface morphologies. For 
this we derive a nonlinear theory that describes the time 
evolution of the surface morphology. At short time scales 
the nonlinear theory predicts the development of a peri- 
odic ripple structure, while at large time scales the sur- 
face morphology may be either rough or dominated by 
new ripples, that are different from those existing at 
short time scales. We find that transitions may take 
place between various surface morphologies as the ex- 
perimental parameters (e.g. angle of incidence, penetra- 
tion depth) are varied. Usually stochastic equations de- 
scribing growth and erosion models are constructed us- 
ing symmetry arguments and conservation laws. In con- 
trast, here we show that for sputter eroded surfaces the 
growth equation can be derived directly from a micro- 
scopic model of the elementary processes taking place in 
the system. A particular case of our theory was pre- 
sented in In addition, we show that the presented 
theory can be extended to describe low temperature rip- 
ple formation as well. We demonstrate that under certain 
conditions ion-sputtering can lead to preferential erosion 
that appears as a surface diffusion term in the equation 
of motion even though no mass transport along the sur- 
face takes place in the system. To distinguish it from 
ordinary surface diffusion, in the following we refer to 
this phenomenon as effective smoothing (ES). We cal- 
culate analytically an effective surface diffusion constant 
accounting for the ES effect, and study its dependence 
on the ion energy, flux, angle of incidence, and pene- 
tration depth. The effect of ES on the morphology of 
ion-sputtered surfaces is summarized in a morphological 
phase diagram, allowing for direct experimental verifica- 
tion of our predictions. A restricted study along these 
lines appeared in |^ . 

The paper is organized as follows. In Section II we 
review the recent advances in the scaling theory of rough 
(self-affine) interfaces. Section III is dedicated to a brief 
overview of the experimental results on surface morphol- 
ogy development under ion sputtering. A short summary 
of the theoretical approaches developed to describe the 
morphology of ion sputtered surfaces is presented in Sec- 
tion IV. This section also contains a short description 
of Sigmund's theory of sputtering, that is the basis for 
our calculations. In Section V we derive the nonlinear 
stochastic equation describing sputter erosion. Analysis 
of this equation is presented in Section VI, discussing 
separately both the high and low temperature ripple for- 
mation. We compare the predictions of our theory with 
experimental results on surface roughening and ripple 
formation in Section VII, followed by Section VIII, that 
summarizes our findings. 



II. SCALING THEORY 

In the last decade we witnessed the development of an 
array of theoretical tools and techniques intended to de- 
scribe and characterize the roughening of nonequilibrium 
surfaces and interfaces |37|. Initiated by advances in un- 
derstanding the statistical mechanics of various nonequi- 
librium systems, it has been observed that the roughness 
of many natural surfaces follows rather simple scaling 
laws, which can be quantified using scaling exponents. 
Since kinetic roughening is a common feature of ion- 
bombarded surfaces, before we discuss the experimen- 
tal results on sputtering, we need to introduce the main 
quantities characterizing the surface morphology. 

Let us consider a two-dimensional surface that is char- 
acterized by the height function h{x,y,t). The morphol- 
ogy and dynamics of a rough surface can be quantified 
by the interface width, defined by the rms fluctuations in 
the height h(x,y,t), 



w{L,t)^J^ J2 [h{x,y)-h]^, 



(1) 



= 1.L 



where L is the linear size of the sample and h is the mean 
surface height of the surface 
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x,y,t). 



(2) 



x,y—l,L 



Instead of measuring the roughness of a surface over 
the whole sample size L x L, we can choose a window of 
size £ X i, and measure the local width, w{£). A general 
property of many rough surfaces is that the roughness 
depends on the length scale of observation. This can be 
quantified by plotting w{i) as a function of £. There are 
two characteristic regimes one can distinguish: 

(i) For length scales smaller than £x, the local width 
increases as 



W{i)r^£", 



(3) 



where a is the roughness exponent. If we are interested 
in surface phenomena that take place at length scales 
shorter than £x then we cannot neglect the roughness of 
the surface. In this regime the roughness is not simply a 
number, but it depends on the length scale accessible to 
the method probing the surface. 

(ii) For £ ^ £x, w{£) is independent of £, thus, at 
length scales larger than £x , the surface is smooth. In this 
regime we can characterize the surface roughness with the 
saturation width Wgat = w{£x)- 

The dynamics of the roughening process can be best 
characterized by the time dependent total width (|l|) . At 
early times the total width increases as w{L,t) ~ t^^, 
where /3 is the growth exponent. However, for finite sys- 
tems, after a crossover time tx, the width saturates, fol- 
lowing the Family- Vicsek scaling function 



wiL,t)^tPg(^±^ 



(4) 
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where z = a//? is the dynamic exponent and g(u 
1, while g{u ^ 1) '--^ . 

ScaUng relations such as Eq. allow us to define 
universality classes. The universality class concept is a 
product of modern statistical mechanics, and encodes the 
fact that there are but a few essential factors that deter- 
mine the exponents characterizing the scaling behavior. 
Thus different systems, which at first sight may appear 
to have no connection between them, behave in a remark- 
ably similar fashion. The values of the exponents a and 
(3 are independent of many "details" of the system, and 
they are uniquely defined for a given universality class. 
In contrast, other quantities, such as A, ix, or Wsat, are 
non-universal, i.e. they depend on almost every detail of 
the system. 



III. EXPERIMENTAL RESULTS 



The morphology of surfaces bombarded by energetic 
ions has long fascinated the experimental community. 
Lately, with the development of high resolution obser- 
vation techniques such as atomic force (AFM) and scan- 
ning tunneling (STM) microscopies, this problem is living 
a new life. The various experimental investigations can 
be classified into two main classes. First, early investi- 
gations, corroborated by numerous recent studies, have 
found that sputter eroded surfaces develop a ripple mor- 
phology with a rather characteristic wavelength of the 
order of a few micrometers P-pTj. However, a number 
of research groups have found no evidence of ripples, but 
observed the development of apparently random, rough 
surfaces p8|~p6|, that were characterized using scaling 
theories. In the following we summarize the key experi- 
mental observations for both ripple development and ki- 
netic roughening. 



A. Ripple formation 

The ripple morphology of ion bombarded surfaces has 
been initially discovered in the mid 1970's |-§. bmce 
then, a number of research groups have provided detailed 
quantitative results regarding the ripple characteristics 
and dynamics of ripple formation. It is beyond the scope 
of this paper to offer a comprehensive review of the vast 
body of the experimental literature on the subject. Thus, 
we selected a few experiments that offer a representative 
picture of the experimental features that appear to be 
common to different materials. 

Angle of incidence: An experimental parameter which 
is rather easy to change in sputtering is the angle of in- 
cidence 9 of the incoming ions relative to the normal to 
the average surface configuration. Consequently, numer- 
ous research groups have investigated the effect of 9 on 
the ripples. These results indicate that ripples appear 



only for a limited range of incidence angles, which, de- 
pending on materials and ions involved, typically vary 
between 30° and 60°. 

Support for a well defined window in 9 for ripple forma- 
tion was offered by Stevie et al. , who observed abrupt 
secondary ion yield changes (correlated with the onset 
of ripple morphology development) in experiments on 6 
and 8 keV sputtering of Si and 8, 5.5, and 2.5 keV 
O2 sputtering of GaAs at incidence angles between 39° 
and 52°. These results were corroborated by Karen et al. 
|^-|lO|, who investigated ripple formation on GaAs sur- 
faces under bombardment with 10.5 keV ions. They 
found that ripple formation takes place for angles of in- 
cidence between 30° and 60° (see Table I of Ref. @). 
Similarly, Wittmaack found that ripple formation oc- 
curs at incidence angles between 32° and 58° during 10 
keV Oj-ion bombardment of a Si target. 

Temperature dependence: Another parameter that has 
been found to influence the ripple structure, and in par- 
ticular the ripple wavelength, is the temperature of the 
substrate, T. Two different behaviors have been docu- 
mented: exponential dependence of the ripple wavelength 
on T has been observed at high temperatures, while the 
wavelength was found to be constant at low tempera- 
tures. 

A series of experiments on the temperature depen- 
dence of ripple formation were reported by MacLaren 
et al. ijl^. They studied InP and GaAs surfaces bom- 
barded with 17.5 KeV Cs+ and 5.5 keV ion beams 
in the temperature range from —50° C to 200° C. For 
GaAs bombarded by Cs+ ions the ripple wavelength in- 
creased from 0.89 /zm to 2.1 /Ltm as the temperature in- 
creased from 0° C to 100° C. Probably the most inter- 
esting finding of their study was that lowering the tem- 
perature, the ripple wavelength did not go continuously 
to zero as one would expect, since the surface diffusion 
const ant d ecreases exponentially with temperature (see 
Sect. tVC), but rather at approximately 60° C it stabi- 
lized at a constant value. MacLaren et al. interpreted 
this as the emergence of radiation enhanced diffusion, 
that gives a constant (temperature independent) contri- 
bution to the diffusion constant. Recently, Umbach et 
al. ||l^ have studied the sputter-induced ripple forma- 
tion on Si02 surfaces using 0.5 — 2.0 keV Ar ion beams. 
The temperature dependence of the ripple wavelength 
£ was investigated for temperatures ranging from room 
temperature to 800° C. It was found that for high tem- 
peratures, T > 400° C, the ripple wavelength follows the 
Arrhenius law (l/T^/^j g^p {-AE/2kBT), indicating the 
thermally activated character of the relaxation processes. 
However, at low temperatures the ripple wavelength was 
independent of temperature, indicating the presence of a 
temperature independent relaxation mechanism. 

Results indicating temperature independent non- 
diffusive relaxation have been reported for crystalline ma- 
terials as well by Carter et al. ||lj]. In these experiments 
Si bombarded with highly energetic 10 — 40 keV Xe"''-ions 
led to ripple formation with wavelength £ ~ 0.4 /im for 
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angles of incidence close to 45°. Changing the surface 
temperature from 100 K to 300 K the ripple wavelength 
and orientation did not change. This observation led 
the authors to conclude that the smoothing mechanism 
is not of thermal origin. They also found that the rip- 
ple wavevector is relatively insensitive to the primary ion 
energy. 

Flux and fluence dependence: The effect of the flux on 
the surface dynamics was studied by Chason et al. ||l5|,|l6| . 
In these experiments a 1 keV Xe ion beam was directed 
towards a Si02 sample with an angle of incidence of 55°. 
The typical incoming flux was 10^^ cm~^s~^ and fluence 
(the number of incoming atoms per surface area, playing 
the role of time) was up to 10 x 10^^ cm^^. The surface 
was analyzed using in situ energy dispersive X-ray reflec- 
tivity and ex situ AFM. It was found that the interface 
roughness, which is proportional to the ripple amplitude, 
increases linearly with the fluence. Similar experiments 
were performed on Ge(OOl) surfaces (l7| using 0.3, 0.5, 
and 1 keV Xe ion beams for T = 350° C. For flux values 
up to 3 /LtA/cm^ and fluences up to 6 x 10^^ cm~^, the 
roughness is seen to increase as the square of the flux. 

Ion energy: The ripple wavelength dependence on the 
incident ion energy and the angle of incidence was re- 
ported in Refs. |p-pi|. The experiments indicate that 
the ripple wavelength is linear in the energy, follow- 
ing £ ~ ecosO. Similar results were obtained in Ref. 
|18| , providing an extensive study of ripple formation by 
secondary ion spectrometry and scanning electron mi- 
croscopy. The ripple topography was observed during 
primary ion bombardment of a Si(lOO) substrate with ion 
energies ranging between 1.5 keV and 9 keV. No ripples 
were observed for energies less than 1.5 keV or for high 
energies, such as 1.5 keV and 7 keV, when Ar"*" primary 
ions were used. The experiments indicate that the ripple 
wavelength increases linearly from 100 to 400 nm when 
the energy of the primary ion changes from 1 to 9 keV. 
Furthermore, the experimental data indicated that the 
primary ion penetration depth a and the ripple wave- 
length £ are related by the empirical relation £ — 40a. 
The wavelength of the ripples is found to be indepen- 
dent of the primary ion flux and dependent merely on 
fluence, i.e. sputtered depth. The recent results by Um- 
bach et al. [ p^ provided further strong evidence for the 
linear relationship between the ion energy and the ripple 
wavelength for Si02 substrates (see below). 

Ripple amplitude: Indirect results on the ripple am- 
plitude were presented by Vajo et al. [p^ in their study 
of the surface topography induced secondary ion yield 
changes on Si02 surfaces bombarded by ions. The 
authors have found that the yield changes exponentially 
in the first stages of ripple development and saturates for 
large sputtered depth. Direct evidence on ripple ampli- 
tude saturation was obtained by Erlebacher et al. po| , 
who measured the time evolution of the ripple amplitude 
in experiments bombarding Si(lOO) surfaces with 0.75 
keV Ar+ ions. They found that, while at short times 
the ripple amplitude increases exponentially, it saturates 



Material 


Ion 
type 


Angle 
or 

incidence 


Ion energy 
[Key ) 


Ripple 
wavelength 
(^m) 


Ref. 


GaAs(lOO) 


oj 


39° 


8 


0.2 






GaAs(lOO) 




42° 


5.5 


0.1 






GaAs(lOO) 




37° 


10.5 


0.23 






GaAs(lOO) 


Ot 


42° 


5.5 


0.13 




ici 

J-V. 1 


GaAs(lOO) 




39° 


8.0 


0.21 





GaAs(lOO) 


Ot 


37° 


10.5 


0.27 




ICl 


GaAs(lOO) 


ot 


57° 


13 


0.33 




GaAs 


ot 


40° 


3.0 


0.075 




2(1 


GaAs 


Ot 


40° 


7.0 


0.130 




mm ' 

1 P ] 
1 


Ge(OOl) 


Xe+ 


55° 


1 


0.2 




Si(OOl) 


Ot 


41° 


6 


0.4 




7 


Si(OOl) 


ot 


42° 


5.5 


0.5 




7 


Si(lOO) 


ot 


39° 


8 


0.5 




7 

ill 


Si(lOO) 


ot 


40° 


3 


0.198 




Si(lOO) 


ot 


40° 


5 


0.302 




1^1 


Si(lOO) 


ot 


40° 


9 


0.408 




li] 

li \ 


Si(lOO) 


Ar+ 


67.5° 


0.75 


0.57 




Si 


Xe+ 


45° 


40 


0.4 




14] 


Si 




37° 


12.5 


0.35 




21] 
1J[ 


SiOa 


Ar+ 


45° 


0.5-2 


0.2-0.55 




SiOa 


Xe+ 


55° 


1 


0.03 




16] 



TABLE I. Summary of the ripple characteristics reported 
in sputter erosion experiments of non-metallic substrates. In 
all cases shown, the ripple wave vector is parallel to the ion 
beam direction. Note that a number of experiments have 
obtained indirect information on ripple formation from sec- 
ondary ion yield changes. These have not been included in 
the table. 

after a crossover time has been reached. Furthermore, 
the experiments indicate that the crossover time scales 
with the temperature induced surface diffusion constant. 

Surface chemistry and other morphological features: 
While a number of attempts have been made to ex- 
plain ripple formation based on chemical effects, such as 
ot variations ||l^,|l],^ , most of these studies were con- 
tradicted by subsequent investigations where such 
chemical component were not present. Furthermore, in 
Refs. [p[-p^ it was unambiguously shown that the process 
of ripple formation is not caused by defects or inherited 
irregularities on the surface, but is determined merely 
by the primary ion characteristics. These results indi- 
cate that ripple formation is independent of microscopic 
details and the surface chemistry. 

Ripple formation on crystalline and metallic surfaces: 
As the discussed experimental results have indicated, rip- 
ple formation takes place under a variety of conditions 
and on surfaces of different materials, including both 
crystalline and amorphous materials. Despite the fact 
that Sigmund's theory, the basis of all theories of ripple 
formation, has been developed for amorphous targets, it 
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is worth noting that these approaches describe many fea- 
tures of ripple formation on crystalline surfaces as well. 
However, when discussing ripple formation on crystalline 
materials, we always have to be aware that additional 
effects, induced by the crystalline anisotropy, could be 
present. An example of ripple development on crystalline 
materials has been obtained for Ag(llO) surfaces under 
low energy (e > 800 eV)Ar+ primary beam bombard- 
ment by Rusponi et al. [M. Ripples with wavelength of 
approximately 600 A, oriented along the (110) crystallo- 
graphic direction, appeared in a temperature range 270° 
K < r < 320° K at normal ion incidence. The ripple 
structure was found to be unstable at room tempera- 
ture, i.e. substantial smoothing of the surface with time 
takes place. The structure depends on the ion flux and 
ion energy. Similar results are available for ion-sputtered 
Cu(llO) monocrystals using a 1 keV Ar+ ion beam [^ . 
For normal incidence a well defined ripple structure was 
observed with wave vectors whose direction changes from 
(001) to (110) as the temperature of the substrate is 
raised. Off-normal sputtering generated ripples whose 
orientation depends both on the ion direction and the 
surface orientation. The authors suggested that this phe- 
nomenon can be explained in terms of anisotropic surface 
diffusion. 

Summary: As the presented results indicate, ripple for- 
mation on ion-sputtered surfaces has been observed by 
many groups in various systems (for a partial summary 
see Table B). The main experimental results, common to 
most studied materials, can be summarized as follows: 

• Off-normal ion bombardment often produces period- 
ically modulated structures (ripples) on the surfaces of 
amorphous and crystalline materials. The wavelength of 
the ripples £ is usually of the order of tenths of microm- 
eters. 

• For non metallic substrates, the orientation of the 
ripples depends on the angle of incidence 9, and in most 
cases is either parallel or perpendicular to the direction 
of the incoming ions. 

• At low temperatures the ripple wavelength is inde- 
pendent of T, while it follows the Arrhenius law £ ^ 
{l/T^/'^)exp{-AE/kBT) at higher temperatures. 

• Numerous experiments find that the ripple wave- 
length is proportional to the ion range, and thus to the 
ion energy for intermediate energies. 

• The ripple wavelength in many cases is independent 
of the ion flux, but systematic flux dependence has also 
been reported. 

• The amplitude of the periodic modulations grows ex- 
ponentially for early times, but saturates after a typical 
crossover time has been reached. In many instances, the 
ripple wavelength £ is found to be independent of time. 

• Evidence for ripple formation was obtained for differ- 
ent materials and different primary ions, suggesting that 
the mechanism responsible for ripple formation is largely 
independent of surface chemistry, chemical reactions on 
the surface, or defects. 



B. Kinetic roughening 



Motivated by the advances in characterizing the mor- 
phology of rough surfaces, recently a number of exper- 
imental studies have focused on the scaling properties 

6l. These 



of surfaces eroded by ion bombardment 
experiments have demonstrated that under certain ion 
bombardment conditions ripples do not form, and the 
surface undergoes kinetic roughening. The goal of the 
present section is to review the pertinent experimental 
results, aiming to summarize the key features that a com- 
prehensive theory should address. 

Surface roughness and dynamical exponents: In the ex- 
periments of Eklund et al. pyrolytic graphite was 
bombarded by 5 keV Ar ions, at an angle of incidence 
of 60°. The experiments were carried out for two ffux 
values, 6.9 x 10^^ and 3.5 x lO"^'* ions s'^ cm~^, the total 
ffuences being 10^^, IQi"^ and 10^* ions • cm"^. STM mi- 
crographs indicated that large scale features develop with 
continuous bombardment, the interface becoming highly 
correlated and rough. The scaling properties have also 
been probed using the Fourier transform of the height- 
height correlation function, obtaining a dynamic expo- 
nent z in the range 1.6 — 1.8, and a roughness exponent 
in the range 0.2 — 0.4. These exponents are consistent 
with the predictions of the continuum theory, describ- 
ing kinetic roughening, proposed by Kardar, Parisi and 
Zhang (KP Z) that predicts z ~ 1.6 and a ~ 0.38 



(see section I V A 1 ) . 

A somewhat larger roughness exponent has been mea- 
sured for samples of iron bombarded with 5 keV Ar 
ions at an angle of incidence of 25° |^0|. The interface 
morphology was observed using STM, and the height- 
height correlation function indicated a roughness expo- 
nent a = 0.53±0.02 The mechanism leading to such 
a roughness exponent is not yet understood in terms of 
continuum theories, since for two dimensional surfaces 
the existing continuum theories predict a values of 0.38, 
2/3 and 1 p7| , far from the observed roughness exponent. 

Anomalous dynamic-scaling behavior of sputtered sur- 
faces was reported by Yang et al. The experi- 
ments performed on Si(lll) surfaces with 0.5 keV Ar+ 
ions with flux 0.2 /lA/mm^ in a wide range of substrate 
temperatures have provided evidence of scaling behav- 
ior in the limit of small distances r. The height-height 
correlation function has been found to follow C(r) = 
{{h{ro)-h{r + ro))^) ^ r^" logt, with a ~ 1.15±0.08for 
temperatures lower than 530° C. No roughening was ob- 
served for higher temperatures, demonstrating the tem- 
perature dependence of kinetic roughening. 

Temperature dependence: The effect of surface relax- 
ation due to surface diffusion on roughening of GaAs(llO) 
surfaces eroded by 2 keV Ar+ and Xe"^ was reported 
by Wang et al. ||3^. They found that both the height- 
height correlation function and the small scale rough- 
ness increase significantly faster during erosion at higher 
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temperatures than at lower ones. The surface width in 
these experiments increased with (3 = 0.3 at T = 725 
K and there was no evidence of scaling for lower tem- 
peratures, such as T = 625 K. The roughness exponent 
has been determined as a = 0.38 ± 0.03. In general, 
Ref. ||3l[] concludes that on large scales the surfaces are 
rougher at higher temperatures, contrary to the expecta- 
tion of smaller roughness due to increased relaxation by 
surface diffusion. Similar conclusions on the temperature 
dependence of the scaling properties were drawn in Ref. 
p3| . A sharp transition between scaling regimes in ion- 
bombardment of Ge(OGl) surfaces with 1 keV Xe ions was 
observed at Tc — 488 K. The regimes above and below Tc 
are characterized by dynamic scaling exponents (3 with 
values 0.4 and 0.1, respectively. The surface roughness 
of Si(lll) during low-energy (500 eV) ion bombardment 
at T = 610 K was studied in Ref. [|4| using STM. It was 
found that the rough morphology is consistent with the 
early time behavior of the noisy Kuramoto-Sivashinsky 
(KS) equation (see Sect. IV A3| ). The measured rough- 
ness exponent was a = 0.7 and the dynamic exponent 
was /3 = 0.25. 

Low energy ion bombardment: Recently a number of 
experiments and simulations have focused on low energy 
ion bombardment (i.e., at energies 50-500 eV), for which 
the secondary ion yields are smaller than one [^-^. 
In this systems, the effect of the ions is limited to the 
surface of the material, the collective effect created by 
the collision cascade being less relevant. Often, such low 
energy sputtering leads to layer-by-layer erosion, almost 
mirroring layer-by-layer growth phenomena. The effect 
of vacancy diffusion and Schwoebel barriers can be rather 
well studied in these systems, that include Ge(OOl) sur- 




, and Si(lll) sur- 
In the absence of 



face etching, by 240 eV Xe ions 
faces etched by 100 eV Ar ions 
the collision cascade, ripple formation and kinetic rough- 
ening seen at higher energies, the subject of this paper, 
do not appear. 

Various experimental results on ion-bombardment in- 
duced surface roughening are summarized in Table ||. 
These experiments demonstrate that kinetic roughening 
is one of the major experimental morphologies generated 
by ion bombardment. However, as Table || indicates, 
there is a considerable scattering in the scaling expo- 
nents. This scattering is not too disturbing at this point: 
accurate determination of the scaling exponents from ex- 
perimental data is rather difficult, since often the scal- 
ing regime is masked by strong crossover effects. As we 
demonstrate later, due to the separation of the linear 
and nonlinear regimes, such crossovers are, indeed, ex- 
pected in sputter erosion. Thus the main conclusion we 
would like to extract from this section is that numerous 
experiments do observe kinetic roughening, and find that 
scaling concepts can successfully characterize the surface 
morphology. It will be a major aim of the theory pro- 
posed here to account for the origin of kinetic roughening 
and predict the scaling exponents. 



Surface 
material 


Ion 

type 


Ion 
energy 
(keV) 


Angle 
of 

incidence 


a 


P 


Ref. 


Graphite 


Ar+ 


5 


60° 


0.2-0.4 


2.5-2.9 




2^1 


Iron 







zo 


U.oo 






3C1 
311 


Si(lll) 


Ar 


0.5 


0° 


1.15* 






Si(lll) 


Ar+ 


0.5 


0° 


0.7 


0.25 




34] 

ii 

3J1 


GaAs(llO) 


Ar+ 


2 


0° 


0.38(3) 


0.3 




Ge(OOl) 


Xe 


1 


30° 




0.1, 0.4 




Ni, Cr, Cu 


Ar+ 


1 


86° 


0.49 






3(1 



TABLE II. Summary of the scaling exponents, character- 
izing the surface morphology, reported in various experiments 
on sputter eroded surfaces. * - anomalous logarithmic scaling 
was reported in this experiment. 



IV. THEORETICAL APPROACHES 

The recent theoretical studies focusing on the char- 
acterization of various surface morphologies and their 
time evolution have revolutionized our understanding of 
growth and erosion (for reviews, see The physi- 

cal understanding of the processes associated with inter- 
face roughening require the use of the modern concepts 
of fra ctal geometry, universality and scaling. In Sect. 
[V A we review the major theoretical contributions to 
this area, necessary to describe th e mor phology of ion- 
eroded surfaces. In Sects. IV B to IVE we then review 



the available theoretical approaches (whether through 
continuum equations or by the use of discrete atomistic 
models) that specifically describe surfaces eroded by ion- 
bombardment, emphasizing the procedures which allow 
to describe within a continuum approach some of the rel- 
evant physical processes taking place at the surface, such 
as surface diffusion and beam fluctuations. 



A. Continuum theories of kinetic roughening 

The full strength of the continuum theories comes from 
the prediction of the asymptotic behavior of the growth 
process valid in the long time and large length scale lim- 
its. These limits are often beyond the experimentally or 
practically interesting time and length scales. A notable 
exception is sputter erosion, where both the short time 
ripple development and the asymptotic kinetic roughen- 
ing have been observed experimentally. Consequently, 
next we discuss separately the continuum theories needed 
to understand sputter erosion. 
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1. Kardar-Parisi-Zhang (KPZ) equation 

The time evolution of a nonequilibrium interface can 
be described by the Kardar-Parisi-Zhang (KPZ) equation 

il 



dh A,„,,n 



(5) 



The first term on the rhs describes the relaxation of the 
interface due to the surface tension {i/ is here a pos- 
itive constant) and the second is a generic nonlinear 
term incorporating lateral growth or erosion. The noise, 
r](x,y,t), reflects the random fluctuations in the growth 
process and is a set of uncorrelated random numbers with 
zero configurational average. For one dimensional in- 
terfaces the scaling exponents of the KPZ equation are 
known exactly, as a = 1/2, (3 = 1/3, and z = 3/2. How- 
ever, for higher dimensions they are known only from 
numerical simulations. For the physically most relevant 
two dimensional surfaces we have a ~ 0.38 and p ~ 0.25 

If A = in Eq. m, the remaining equation describes 
the equilibrium fluctuations of an interface which tries 
to minimize its area under the influence of the external 
noise. This equation, first introduced by Edwards and 
Wilkinson (EW) [H, can be solved exactly due to its lin- 
ear character, giving the scaling exponents a = (2 — o?)/2 
and P = (2 — c?)/4. For two dimensional interfaces {d — 2) 
we have a — f3 — 0, meaning logarithmic roughening of 
the interface, i.e., 'w{L) ~ logL for saturated interfaces, 
and w(t) ~ \ogt for early times. 



2. Anisotropic KPZ equation 

The presence of anisotropy along the substrate may 
drastically change the scaling properties of the KPZ 
equation. As a physical example consider an ion bom- 
barded surface, where the ions arrive under oblique in- 
cidence in the x — h plane. As a result, the x and y di- 
rections along the substrate will not be equivalent. This 
anisotropy is expected to appear in the erosion equation, 
leading to an anisotropic equation of the form {d = 2) 



dh 



= J.^dih + Vyd'h+^{d^h) 



dt 



+^{dyhf +r]{x,y,t), 



[AKPZ] 



(6) 



where dxh = dh/dx and dyh = dh/dy. The anisotropy 
leads to surface tension and nonlinear terms that are 
different in the two directions, which have been incor- 
porated in the growth equation by considering different 
values for the coefficients v and A (in Eq. (^), i>x and 
Vy are positive constants). Equation (|^) is called the 
anisotropic KPZ (AKPZ) equation. It was introduced 
by Villain 152] , and its nontrivial properties were studied 



by Wolf ]|5|,Q. We note that ii = v_u and A^; = Aj^, 
Eq. d) reduces to the KPZ equation (|). The AKPZ 
equation has different scaling properties depending on 
the signs of the coefficients A^; and Ay. When A^; • Aj, < 0, 
a surface described by the AKPZ equation has the same 
scaling properties as the EW equation. However, when 
As • Ay > the scaling properties are described by the 
isotropic KPZ equation (||). Thus, changing the sign of 
Aa; or \y can induce morphological phase transitions from 
power law scaling {w ~ t^;w{L) L") to logarithmic 
scaling {w ^ logt; w{L) ~ logL). 



3. Kuramoto-Sivashinsky (KS) equation 

The Kuramoto-Sivashinsky (KS) equation, originally 
proposed to describe chemical waves and flame fronts 
p5[ , is a deterministic equation of the form: 



WW^h- 



X 



K\7^h+ -{Vhy [KS]. 



(7) 



dh 

While it is deterministic, its unstable and highly nonlin- 
ear character gives rise to chaotic solutions. The analysis 
of the KS equation for one dimensional surfaces shows 
]^6|-|62| that in the limit of long time and length scales, 
the surface described by the KS equation is similar to 
that described by the KPZ equation, i.e. obeys self-affine 
scaling with exponents z = 3/2 and (3 = 1/3. The short 
time scale solution of KS equation reveals an unstable 
pattern-forming behavior, with a morphology reminis- 
cent of ripples Q. For two dimensional surfaces, how- 
ever, the results are not clear. Computer simulations are 
somewhat contradictory, providing evidence for both EW 
and KPZ scahng 

The anisotropic KS equation was studied in Ref. [ |65[ , 
indicating that for some parameter values the nonlin- 
earities cancel each other, and lead to unstable modes 
dominating the asymptotic morphology. At early times 
the surface displays a chaotic pattern, with stable do- 
mains that nucleate and grow linearly in time until ripple 
domains of two different orientations are formed. The 
pattern of domains of perpendicularly oriented ripples 
coarsen with time until one orientation takes over the 
system. 

There are various physical systems, including ion sput- 
tering, in which the relevant equation for the surfac e 
height is a noisy version of the KS equation 
Dynamical renormalization group analysis [ 
surface dimensions d = I and 2 indicate that the large 
distance and long time behavior of such noisy generaliza- 
tion of Eq. (0) is the same as that of the KPZ equation, 
the d = 2 result being only quantitative. 



B. Bradley and Harper theory of ripple formation 

A rather successful theoretical model, capturing many 
features of ripple formation, was developed by Bradley 



©il- 

3§ for the 
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and Harper (BH) p8[ . The y used Sigmund's theory of 
sputtering (see Sect. IV E ) to relate the sputter 

yield to the energy deposited onto the surface by the in- 
coming ions. This work has demonstrated for the first 
time that the yield variation with the local surface cur- 
vature induces an instability, which leads to the forma- 
tion of periodically modulated structures. This instabil- 
ity is caused by the different erosion rates for troughs 
and crests, the former being eroded faster than the lat- 
ter (see Fig. 0). Consequently, any surface perturbation 
increases exponentially in time. Viewing the surface pro- 
file as a smooth analytical function of coordinates, BH 
assumed that the surface can be locally approximated by 
a quadratic form. Due to the erosion mechanism, de- 
scribed in Fig. |l|, the erosion rate depends on the local 
curvature. Combining the curvature dependent erosion 
velocity with the surface smoothing mechanism due to 
surface diffusion (see and next section), BH de- 

rived a linear equation for surface morphology evolution 




a) convex 




b) concave 



vy{e)dlh- KV^h. (8) 



Here Vx{0), Vyid) are the effective surface tensions gener- 
ated by the erosion process, dependent on the angle of 
incidence of the ions, 9, K is the relaxation rate due to 

surface diffusion [K = Dsjfl^n/kBT exjp | | , where 

AE is the activation energy for surface diffusion, 7 is the 
surface free energy per unit area, T is temperature, Dg 
is the surface diffusion constant, is the atomic volume 
and n is the number of molecules per unit area on the 
surface). The physical instability illustrated in Fig. 1 
leads to the negative signs of the , i^y coefficients in Eq. 
(^. Eq. (H) is linearly unstable, with a Fourier mode kc 
whose amplitude exponentially dominates all the others. 
This mode is observed as the periodic ripple structure. 
Using linear stability analysis, BH derived from Eq. (||) 
the ripple wavelength as 




(JT) 



-1/2 



exp 



(9) 



FIG. 1. Schematic illustration of the physical origin of the 
instability during ion erosion of nonplanar surfaces. A surface 
element with convex geometry (a) is eroded faster than that 
with a concave geometry (b), due to the smaller distances 
(solid lines) the energy has to travel from the impact point to 
the surface (A or A' points). 

Furthermore, the BH model cannot account for low 
temperature ripple formation since the only smoothing 
mechanism it considers is of thermal origin. At low tem- 
peratures the ion energy and flux dependence of the rip- 
ple wavelength also disagree with the BH predictions. 
Despite these shortcomings, the BH theory represents a 
major step in understanding the mechanism of surface 
evolution in ion sputtering since for the first time it un- 
covered the origin of the ion induced surface instability. 
Recently a generalization of BH linear theory has been 
successfully introduced |Q to account for the thermally 
activated anisotropic surface diffusion present in metallic 
substrates such as Cu(llO). 



C. Surface diffusion and deposition noise 



where v is the largest in absolute value of the two neg- 
ative surface tension coefficients, and fj,, and J is 
the ion flux. The calculation also predicts that the rip- 
ple direction is a function of the angle of incidence: for 
small 9 the ripples are parallel to the ion direction, while 
for large 9 they are perpendicular to it. As subsequent 
experiments have demonstrated ||l^,|l^, the BH model 
predicts well the ripple wavelength and orientation. On 
the other hand, the BH equation (^) is linear, predicting 
unbounded exponential growth of the ripple amplitude, 
thus it cannot account for the stabilization of the rip- 
ples and for kinetic roughening, both phenomena being 
strongly supported by experiments (see Sect. Ill A 111 B). 



At high temperatures surface diffusion and fluctuations 
in the ion beam flux are relevant physical mechanisms 
taking place on the surface |Q. In this section, we dis- 
cuss the standard approach to include these phenomena 
in continuum models. Let us consider the simplest sce- 
nario: atoms are deposited on a surface, whereupon they 
diffuse. If we assume that surface diffusion is the only 
relaxation mechanism present, the height h obeys a con- 
tinuity equation of the form 

dh ~ 

^+V.,=0, (10) 

where j is a surface current density tangent to the sur- 
face, and V is calculated in a frame with axes parallel to 
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the surface |7^. In general, j is given by the gradient of 
a chemical potential ^, 



(11) 



where /i minimizes the free energy functional of the sur- 
face !F[h] and is the surface Laplacian or the Laplace- 
Beltrami operator. Taking the latter to be proportional 
to the total surface area 



(12) 



with g as defined in Appendix and neglecting third or 
higher powers of derivatives of h one arrives |73| at 



dt ^ ^ 



(13) 



where K is a. positive constant. Eq. (|l3| ) is the so- 
called linear MBE equation |Q . For an amorphous solid 
in equilibrium with its vapor Eq. was obtained in 
1^,^, together with the expression for the coefhcient 
K as in Eq. §). 

In addition to the deterministic processes, there is con- 
siderable randomness in sputter erosion due to fluctu- 
ations in the intensity of the ion beam. The ion flux 
is defined as the number of particles arriving on the 
unit surface (or per lattice site) in unit time. At large 
length scales the beam flux is homogeneous with an aver- 
age intensity J, but there are local random fluctuations, 
r]{K,t) = SJ{x,t), uncorrelated in space and time. We 
can include fluctuations in Eq. (|l^ ) by considering the 
ion flux to be the sum of the average flux J and the noise 
77, which has zero average, 



(77(x,i)) =0 



and is uncorrelated. 



(77(x,t),7(x',i'))- J'5(x-x')(5(t-t'), 



(14) 



(15) 



where we have assumed a Poisson distribution for the 
shot noise. Consequently, the stochastic growth equa- 
tion describing surface diffusion and fluctuations in an 
erosion process has the form 



dh 
'dt 



(16) 



This variant of Eq. (O) was introduced independently 
by Wolf and Villain and by Das Sarma and Tam- 
borenea fn^ , and played a leading role in developing our 
understanding of MBE. We will use the methods leading 
to (^6|) to incorporate the smoothing by surface diffusion 
in our model of ion erosion. Note, however, that as nu- 
merous experimental studies ||7^-|8^ indicate, ion bom- 
bardment leads to an enhancement of the surface adatom 
mobility and thus may drastically change the relaxation 
mechanism, as compared to regular surface diffusion. 



D. Microscopic models of ripple formation and 
roughening 



Computer simulations provide invaluable insight into 
microscopic processes taking place in physical systems. 
Consequently, a number of recent studies have focused 
on modeling ripple formation at the microscopic level. 
These studies have proven useful in resolving issues re- 
lated to low temperature ripple formation and provided 
important ideas regarding the physical mechanism gov- 
erning ripple formation [p4| p3[ . Here we shortly discuss 
the conclusions reached in some of the most representa- 
tive numerical work. 

Monte Carlo simulations of sputter-induced roughen- 
ing were reported by Koponen et al. |]85|-p0t . Roughening 
of amorphous carbon surfaces bombarded by 5 keV Ar"*" 
ions was studied in [^5|-p7| for incidence angles between 
0° and 60°. It was found that ion bombardment induces 
self-affine topography on the submicrometer scale, the 
roughness exponent being a ~ 0.25 — 0.47 depending on 
the angle of incidence ||86| , |87|] . The growth exponent j3 
was found to be strongly dependent on the relaxation 
mechanism used and changed from /3 ~ 0.3 in the model 
without relaxation to /3 ~ 0.2 — 0.14 when different re- 
laxation rules were used in the simulations. At the same 
time the roughness exponent a was found to be relatively 
insensitive to the relaxation process on the nanometer 
scales. Analogous results were obtained for C ions 
In this Reference, the ripple wavelength was found to be 
relatively independent of the ion energy or the magnitude 
of surface diffusion. Ripple formation was observed even 
at zero temperature, when surface diffusion was switched 
off, indicating the presence of ion induced smoothing. 
Furthermore, these simulations led to the observation of 
travelin g rippl es, as predicted by continuum theories (see 
section VI A 1 ) . Similarly, for 5 keV Ar+ bombardment 
of amorphous carbon substrates, the ripple wave vector is 
seen to change from parallel to normal to the beam 
direction as the incidence angle is increased, in agree- 
ment with BH linear theory, (see Sect. IV). The ripple 
structure was again observed even when no explicit re- 
laxation mechanism was incorporated in the simulations, 
and ripple travelling also occurs. For length scales com- 
parable to the cascade dimensions, self-affine topography 
is observed. 

A discrete stochastic model was introduced in Ref. 
[ pl| , p2t to study the morphological evolution of amor- 
phous one dimensional surfaces under ion-bombardment. 
This is a solid-on-solid model incorporating the erosion 
rate dependence of surface curvature, the local slope de- 
pendence of the sputtering yield, and thermally activated 
surface diffusion. Up to four different dynamical regimes 
have been identified. Initially the surface relaxes by sur- 
face diffusion with a growth exponent (3 ~ 0.38, until 
the onset of the linear BH instability. The instability in- 
duces rapid growth {(3 > 0.5). In this regime the local 
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slopes increase rapidly, which triggers non-linear effects 
eventually stabilizing the surface, (3 taking up the EW 
value /3 ~ 1/4, which indicates that an effective positive 
value of the surface tension has been generated. Finally, 
in the asymptotic time limit [3 reaches the KPZ value 
/3 ~ 1/3. This behavior is consistent with that displayed 
by the noisy KS equation [Q. Furthermore, the ana- 
lytical study 1 94 1 using the master equation approach to 
interface models has shown that the noisy KS equa- 
tion indeed provides the continuum limit of the discrete 
stochastic model of Ref . . Conversely, the results of 
the simulations in support the theoretical conclusions 
of Ref. Q that the asymptotic behavior of the noisy KS 
equation is the same as that of the KPZ equation for one 
and two dimensional surfaces. 

In summary, Monte Carlo simulations of the sputtering 
process of amorphous materials have shown that interme- 
diate and high energy ion bombardment may lead to rip- 
ple formation in a wide parameter range. Furthermore, 
ripple formation was observed even at zero temperature. 
Computer simulations have also confirmed the linear de- 
pendence of the ripple wavelength on the incident ion 
penetration depth and the fact that ripple formation is 
a process fully determined by the incident ion charac- 
teristics and not caused by any defects, irregularities or 
surface chemistry. The same simulations have confirmed 
that under some bombarding conditions the surface is 
rough, and obeys scaling. 



E. Sigmund's theory of sputtering 

The erosion rate of ion bombarded surfaces is charac- 
terized by the sputtering yield, F, defined as the average 
number of atoms leaving the surface of a solid per in- 
cident particle. In order to calculate the yield and to 
predict the surface morphology generated by ion bom- 
bardment, we first need to understand the mechanism of 
sputtering, resulting from the interaction of the incident 
ions and the substrate 0J^. In the process of sputter- 
ing the incoming ions penetrate the surface and transfer 
their kinetic energy to the atoms of the substrate by in- 
ducing cascades of collisions among the substrate atoms, 
or through other processes such as electronic excitations. 
Whereas most of the sputtered atoms are located at the 
surface, the scattering events that might lead to sputter- 
ing take place within a certain layer of average depth a, 
which is the average penetration depth of the incident 
ion. A qualitative picture of the sputtering process is 
as follows; an incoming ion penetrates into the bulk of 
the material and undergoes a series of collisions with the 
atoms of the substrate. Some of the atoms undergo sec- 
ondary collisions, thereby generating another generation 
of recoiling atoms. A vast majority of atoms will not 
gain enough energy to leave their lattice positions per- 
manently. However, some of them will be permanently 
removed from their sites. 




FIG. 2. Schematic illustration of the energy distributed by 
an incident ion. While the collision process induced by a ion 
is rather complex, according to Sigmund it can be reduced to 
the following effective process: The ion penetrates the bulk of 
the material and stops at point P, where all its kinetic energy 
is released and spread out to the neighboring sites following 
a Gaussian form with widths o and ^l. 

The atoms located in the close vicinity of the surface, 
which can gain enough energy to break their bonds, will 
be sputtered. Usually the number of sputtered atoms 
is orders of magnitude smaller than the total number of 
atoms participating in the collision cascade. 

A quantitative description of the process of ion sputter- 
ing was developed by Sigmund | |69[ |. Assuming an amor- 
phous target, Sigmund derived a set of transport equa- 
tions describing the energy transfer during the sputter- 
ing process. A practically important result of Sigmund's 
theory is the prediction of the deposited energy distribu- 
tion: the ion deposited at a point P inside the bulk of 
the material spreads its kinetic energy according to the 
Gaussian distribution 

20-2 



E{r') 



(2^)3/VAi' 



■ exp 



X 



12 



Y 



n 



2p? 



(17) 



In (|T^) Z' is the distance from point P to point O mea- 
sured along the ion trajectory, and X\ Y' are measured 
in the plane perpendicular to it (see Fig. ^ and the inset 
of Fig. H); e denotes the total energy carried by the ion 
and a and /i are the widths of the distribution in direc- 
tions parallel and perpendicular to the incoming beam 
respectively. Deviations of the energy distribution from 
Gaussian ( |l7| ) occur mainly when Mi > M2 , where Mi is 
the mass of the projectile and M2 is the mass of the tar- 
get atom. As shown by Sigmund and Winterbom 
Pql , electronic stopping doesn't affect much the shape 
of deposited-energy distribution. Subsequently, Monte 
Carlo simulations of the sputtering process have demon- 
strated that the deposited-energy distribution and dam- 
age distribution can be well approximated by Gaussian 
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for intermediate and high energies. In general, compar- 
ison of Sigmund's theory with experimental results has 
shown that it describes well the qualitative behavior of 
the backsputtering yield, and in many cases there is good 
quantitative agreement as well ||. 

A quantity of central importance is the mean path 
length of an ion traveling inside the bulk of the mate- 
rial (see Fig. ||), often called penetration depth, given 
by 



aie) = 



1 — TO 



2 TO 



-7 



^2m 



NC„ 



(18) 



where N is the target atom density, Cm is a constant 
dependent on the parameters of the interatomic poten- 
tial JtoI and the exponent to = m{e) varies slowly from 
TO = 1 at high energies to to ~ at very low energies. 
In the region of intermediate energies, i.e. for e between 
10 and 100 keV, to ~ 1/2 and we can approximate the 
penetration depth as a(e) ~ e. 

Eq. ( p^ ) describes the effect of a single incident ion. 
Actually, the sample is subject to an uniform flux J of 
bombarding ions, penetrating the solid at different points 
simultaneously, such that the erosion velocity at an arbi- 
trary point O depends on the total power So contributed 
by all the ions deposited within the range of the distribu- 
tion (^7|) . If we ignore shadowing effects and redeposition 
of the eroded material, the normal erosion velocity at O 
is given by 



Vo =p / dr $(r) E{r), 
Jtz 



(19) 



where the integral is taken over the region TZ of all points 
at which the deposited energy contributes to £o, ^{r) is 
a local correction to the uniform flux J due to variation of 
the local slopes, and the material constant p depends on 
the surface binding energy and scattering cross-section 
MM as 



1 



P 



47r2 NUoCo ' 



(20) 



where Uo is the surface binding energy and Co is a con- 
stant proportional to the square of the effective radius of 
the interatomic interaction potential. 

While the predictions of Sigmund's theory have been 
checked on many occasions, it also has well known lim- 
itations. Next we list two, that will limit our theory 
on the surface morphology as well: (a) the assumption 
of random slowing down and arbitrary collisions works 
satisfactorily only at intermediate and high energies, i.e. 
when e ~ 1 — 100 keV, but may break down at low ener- 
gies; (b) the assumption of a planar surface may influence 
the magnitude of the yield, since surface roughness has 
a tendency to increase the yield 1 97 9^] . 



V. CONTINUUM EQUATION FOR THE 
SURFACE HEIGHT 



Sigmund's theory, while offering a detailed description 
of ion bombardment, is not able to provide direct in- 
formation about the morphology of ion-sputtered sur- 
faces. While Eq. (^9|) provides the erosion velocity, in 
the present form it cannot be used to make analytical 
predictions regarding the dynamical properties of sur- 
face evolution. To achieve such a predictive power, we 
have to eliminate the nonlocality contained in the in- 
tegral (|l9|) and derive a continuum equation describing 
the surface evolution depending only on the local surface 
morphology. The main goal of this section is to provide 
a detailed derivation of such an equation starting from 
Eq. ([l^) . The properties of the obtained equation will be 
discussed in the following sections. 

We start by summarizing the main steps that we follow 
in the derivation of the equation for the surface morphol- 
ogy evolution: 



(i) Using Eq. (19), we calculate the normal component 
of the erosion velocity Vq at a generic point O of the sur- 
face. This calculation can be performed in a local frame 
of reference {X,Y,Z), defined as follows: the Z axis is 
chosen to be parallel to the local normal to the surface 
at point O. The Z axis forms a plane with the trajec- 
tory of an ion penetrating the surface at O. We choose 
the X axis to lie in that plane and be perpendicular to 
Z. Finally, Y is perpendicular to the {X, Z) plane and 
completes the local reference frame, as shown in Fig. |[ 
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FIG. 3. Illustration the local reference frame {X,Y,Z). 
The Z axis is parallel to the local normal to the surface h. 
The ions arrive to the surface along — m. The X axis is in the 
plane defined by Z and m, while the Y axis is perpendicular 
to this plane. The laboratory coordinate frame {x, y, h) has 
its h axis perpendicular to the flat substrate, h and m define 
the {x, h) plane and y is perpendicular to it. The incidence 
angle measured in the local reference frame is 99, and 6 in the 
laboratory frame. 
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(ii) Next we relate the quantities measured in the lo- 
cal frame {X, Y, Z) to those measured in the laboratory 
frame (i, y, h). The latter is defined by the experimental 
configuration as follows: h is the direction normal to the 
uneroded flat surface. The ion direction together with 
the h axis define the (x, h) plane. Finally, the y axis 
is perpendicular to the [x, h) plane (see Fig. || and Ap- 
pendix Furthermore, we have to take into account 
the fact that the local angle of incidence tp, which is the 
angle between the ion trajectory and the local normal to 
the surface, changes from point to point along the sur- 
face. Consequently, tp \s a, function of the local value of 
the slope at O (as measured in the laboratory frame) , and 
the angle 9 between the ion trajectory and the normal h 
to the uneroded surface. 

(iii) Finally, to obtain the equation of motion for the 
surface profile h{x,y,t), we have to project the normal 
component of the velocity of erosion onto the global h 
axis. The time derivative of h{x,y,t) at any point O on 
the surface is proportional to the surface erosion velocity 
Vo at that point and the local normal is defined by the 
gradient of the surface profile h{x, y, t) at O. 

Having defined our objectives and outlined the strat- 
egy, we move on to the description of the calculations. 
We consider point O to be the origin of the local system 
of coordinates {X, Y, Z). To describe the surface profile 
in a neighborhood of O we assume that the surface can 
be described in terms of a smooth analytical infinitely 
differentiable function, i.e. there are no singularities and 
overhangs, and thus we can approximate the surface pro- 
file at an arbitrary point {X, Y, Z) by |93| 



Z{X,Y) 



E 



A, 



-jm+n— 1 



(21) 



m,n— 0,n+m— 3,4 



where, for later convenience, we introduced the following 
notations: 



jn+m-l d'^+rn^f^x, Y) 



1 'm ' 



dnXd'^Y 



(22) 



Here A20 and A02 are proportional to the principal cur- 
vatures of the surface, i.e., to the inverses of the principal 
radii of curvature, Rx and Ry- It must be noted that, 
in our approximation, X and Y (see Fig. ^) are the two 
principal directions of the surface at O, along which the 
curvatures are extremal. This implies the absence in Eq. 
( pT| ) of cross-terms of the type ^ XY , i.e., we neglected 
the term Z{X,Y) / dXdY at O. 

Due to its exponential nature, the deposited energy 
distribution (p7| ) decays very fast and, consequently, 
only particles striking the surface at a point {X, Y, Z) 
such that X/a, Y/a are of order unity, contribute non- 
negligibly to the energy reaching the surface at O. We 
further assume that the surface varies slowly enough so 



that Rx , Ry and the inverses of the higher order deriva- 
tives are much larger than the penetration depth a, i.e. 
the surface is smooth on length scales close to a (this fact 
is supported by nearly all experimental results). Now we 
can calculate the various factors appearing in the integral 

To proceed with Eq. ( |19| ) we note that, with respect 
to local surface orientation, only the normal component 
of the incident flux contributes to ion erosion. Figure 
^ illustrates the calculation of the normal component of 
the flux. In the flgure we consider a point at the surface 
{X,Y,Z), and two other points also on the surface, at 
inflnitesimal distances L and N away from the former. 
We can estimate the correction to the average flux J due 
to the surface slopes by projecting a square perpendicu- 
lar to the ion beam with area n x / onto the surface area 
element intersected by the ion trajectories. The result is 



$(r) ~ J 



/ n 



where J is the average flux. From Fig. ^ 

'dZ\ dZ 



Af — tan ^ 



dX 



dX' 



and 



^ / A ^ 9Z . 

— = COS((^ -I- A(y9j ~ COS (fi — -J— Sm If. 

L uX 
On the other hand, we also have (see Fig. |^) 



n 
N 



cos a 



1 - 



1 fdZ 



2\dY 



1, 



(23) 



(24) 



(25) 



(26) 



so that, combining Eqs. (p3|)-(Pq), and neglecting powers 
of derivatives of the height, we obtain the correction to 
the flux 



<i>(r) ~ J(cos(p -f {dxZ) sinv?). 



(27) 




FIG. 4. Illustration of the calculation of the local correc- 
tion to the average flux J due to the surface curvature. 
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where we used the following notations 



FIG. 5. Reference frame for the calculation of the erosion 
velocity at point O. Following a straight trajectory (thick 
solid line) the ion penetrates an average distance a inside the 
solid (dotted line) after which it completely spreads out its 
kinetic energy. The energy released at point P contributes to 
erosion velocity at O. Inset shows the lateral view for Y' — 0. 





al siiKfi, 


(31a) 


Bi = 


al sin^ ip + a'^ cos^ (p, 


(31b) 


B2 = 


2 

a^ cos 


(31c) 




- al)sm(pcos(p, 


(31d) 


D ^ 




(31e) 


L = 


V2- 


(31f) 



It must be noted that Eq. (|3C|) coincides with the two- 
dimensional version of the local erosion velocity derived 
in Ref. Now we use the approximation for the sur- 
face profile given by Eq. (pl|). Taking a small A„m (see 
Eq. (22)) expansion of the C and B2 coefficients in Eq. 
( pO| ) and evaluating the Gaussian integrals over (^x and 
Cy, we get 



Vo = exp {~al/2) exp — 



1 



X [cos ip + A20 + ro2 A02 + A30 
+ r2lA2ir4oA4o + r22A22 + ro4Ao4] . 



(32) 



Within the same approximation, the surface element 
dr in Eq. ( [i9| ) can be obtained in the local coordinate 
system {X,Y,Z) as 

dr ~ dXdY. (28) 

Next we determine the distances X' , Y' , Z' appearing 
in the exponential distribution (|l^), evaluating them in 
the local coordinate system. Using Fig. ^, we have 

X' — X cos ip + Z sin ip, 
Y' = Y, 

Z' — a + X sin ip — Z cos ip. 

(29) 

Using these expressions, the correction to the ion flux 
(p7|), the deposited energy distribution (|l^) and the ex- 
pression for the surface area element dr (^8|), we can 
calculate the integral ( p9| ) providing the erosion velocity 
Vo- Introducing the dimensionless variables Cx = X/a, 
(■y = Y/a^ and Cz = Z/a, and extending the integration 
limits to infinity, we obtain the following expression for 
the erosion velocity in the laboratory coordinate frame 



Vo 



■exp(-a2/2) 



CT/^2(27r)3/2 

dC,x dC,Y (cos ip 



dCz 



X 



sin ip) 



X exp(-Cyi 



lexp(-CxA)exp(--BiCi) 



X exp(-4i?Cl)cxp(-2CCxCz)exp(B2Cz), (30) 



The expressions for the coefficients Tnm can be found in 
Appendix ^ 

Next we have to rewrite Vo in terms of the labora- 
tory coordinates {x, y, h), which we perform in two steps. 
First, we write the angle ip as a function of 9 and the 
slopes of the surface at O as measured in the laboratory 
frame. Second, we perform the transformation between 
the local and the laboratory coordinates. For both steps 
we will have to make expansions in powers of deriva- 
tives of h{x,y,t). In line with our earlier assumption 
on the smoothness of the surface, we will assume that h 
varies smoothly enough so that we can neglect products 
of derivatives of h for third and higher orders. In the 
laboratory frame, the neglect of overhangs allows us to 
describe a generic point at the surface, such as O, with 
coordinates (x, y, h{x, y)). Considering now the unit vec- 
tors n, m shown in Fig. |^, the angle p is given by 



cos ip = m ■ n 



cos 9 — [dxh) sin 6* 



sin ip = (sin^ 9 + 2{dxh) sin 9 cos 9{dxhf cos^ 9 



2\-l/2 



(33a) 



(33b) 



Thus far, expressions (33a)-(33b) are exact, and the val 



ues of dxh and dyh are already evaluated in the labo- 
ratory frame of reference. To implement our approxi- 
mations, in principle we have to separate the cases for 
normal [9 = 0) and off-normal {9 ^ 0) incidence. Nev- 
ertheless, it can be shown that the former case can be 
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obtained as a smooth limit of the latter. Therefore in 
the following we give the expressions pertaining to the 
ofF-normal incidence and refer the reader to A ppen dix O 
for details on the 6 — limit. Expanding (33a) and ( ^3b| ) 
in powers of the surface height derivatives, we obtain 



cos (fi ~ cos 6 — {dxh) sin 6 

-^{{d^hf + {dyhf) cos 9, 

sin 93 ~ sin 6* + {dxh) cos 6* — ^{dxh)^ sin 9 
1 ,„ , .tCOs^ 9 



(34a) 



(34b) 



Note that these expressions are invariant under the coor- 
dinate transformation y — > — y, but not under x — *■ —x, 
a consequence of 9 being non-zero and of our choice of 
coordinates. Naturally, the 9^0 limit restores the sym- 
metry in the x direction. 



Having obtained the expressions (34a) and (34b), we 
can return to Eq. ( ^2| ) to calculate the dependence of Vq 
on the slopes at O. Finally, all derivatives in (^2|) have 
to be expressed in terms of the laboratory coordinates. 
This can be accomplished given the relation between the 
base vectors of the local frame (X, Y, Z) and those of the 
laboratory frame (x, y, h), derived in Appendix If the 
coordinates of a generic vector r are given by 

r = XX + YY + ZZ local frame, 

r ~ XX + yy + zz laboratory frame, (35) 

then these quantities are related to each other through 



= M 




(36) 



where M is a matrix which has as columns the com- 
ponents of the (X,Y,Z) set of vectors in terms of the 
{x,y,z) (see Appendix^. To obtain the expression for 



J 



dh 

dt 



dh 

OX 

d^h 



dx^dy^ 



dh\ fd'^hX 
dx J \ dx"^ J 

d^h 



D„ 



d^h 



the erosion velocity, ana log ous to ( |32| ) , in the laboratory 
frame, we use Eqs. (34a), ( ^4b| ), and M. along with the 
chain rule for differentiation, and perform expansions in 
powers of derivatives of h{x^ y, t). After some algebra we 
obtain in the laboratory frame 



(37) 



up to fourth order in products of derivatives of h{x,y). 
To summarize our results thus far, within the approxi- 
mations leading to (^2|) and neglecting nonlinearities of 
cubic and higher orders in derivatives of h in the labora- 
tory frame, we obtained Eq. (|37|), providing the relation 
between t he d erivati ves in the two reference frames, and 
relations (34a) and (34b) for the angle of incidence mea- 
sured in the local frame as a function of the angle of 
incidence 9 and of the surface slopes. 

Finally, to relate the velocity of erosion Vq, which is 
normal to the surface at O, to the velocity of erosion of 
the surface along the h axis, dh/dt, we have to project 
the former onto the latter, obtaining 



dh 
'dt 



(38) 



where the negative sign accounts for the fact that Vq 
is the rate at which the surface is eroded, i.e. the aver- 
age height decreases. Furthermore, taking into account 
surface diffusion effects, together with the fluctuations 
(shot noise) in the fl ux of the bombarding particles, as 

we complete 



discussed in Section IV C 
these physical effects 



^ by adding 



dh 
'dt 



-VoV9-KV^h + 7j{r,t), 



(39) 



Finally, wc have to write down the contribution of the 
— Vby^ term to the evolution equation ( p9|) . P erfo rming 
a small slope expansion and using Eqs. (|32|), (34a), and 
(p4q), we obtain 



dh 
dx 

KV^h 



d^h\ 


d^h 


+ Vy 


d^h 


dy^J 




dy^ 




fdhV 


\y 


fdh 


+ T 




2 


\dy 



d^h 
dx'^ 



d^h 
dxdy"^ 



(40) 



where the coefffcients are given by the expressions 



Wo = Fc, 



(41) 



2 „2 2 



(42) 



Fa 



2/3 ^^""^ 



4222, 2222 

a„a„s c +a„a„s c 



4 41 



(43) 



14 



(44) 



-Fa- 



2/ 



(45) 
(46) 



aisc 



atapc'sH^ + 3.^) - aUpcHl + s') + ala^i^ - is') - ia^) , 



(47) 



4 = Fa^-0 {-aUy + ay {2 + c') - ay + ayy{3 - 2.^)} , 



(48) 



(49) 



0,2 = Fa- 



1 1 

6/4 



{-isfif + ay) + ay{3alsf + ay)f + 2{al - al)c\Zfs + ^ay + ay)) , 



(50) 



1 



y^iiayf + ay)f + ay {if + eayj + ay)f 



24 f 

+ 2{al - al)c'il5ayf + lOayf + ai^y , 



(51) 



24 f al ' 



(52) 



= 7^f| y^i<^>')f + <^wy + ayf) + 2(ai ayyayj + ay)} 



(53) 



r 



In the above expressions, we have defined 
F 



f7/i\/27r/ 

and, as introduced in Appendix pi 



_ a _ a 
a^j — , flu — , 
cr /i 



s = smfy, c = cos( 

J. _ 2 2 I 2 2 

f = a„s + a^c . 



(54) 



(55) 



Equation (gg) with the coefficients (^)-(|5|), fuUy de- 
scribe the nonhnear time evolution of sputter eroded sur- 
faces, provided that the leading relaxation mechanisms 
are thermally activated surface diffusion and ion-induced 



effective smoothing. Due to its highly nonlinear char- 
acter, Eq. ( pO[ ) can predict rather complex morphologies 
and dynamical behaviors. In the remainder of the paper 
we will focus on the physical interpretation of the co- 
efficients (|4l|)-(|5^), uncovering their dependence on the 
experimental parameters, and we discuss the morpholo- 
gies described by Eq. (|4(]|). Consistent with the sym- 
metries imposed by the geometry of the problem, the 
coefficients in Eq. ( 40 ) are symmetric under the transfor- 
mation y —y but not under x —x, while for 6 
the system is isotropic in the x and y directions, specifi- 
cally ^ = ^ S,y = ^1 ^ ^2 = 0, Xx = Xy, Vx = Vy, and 

F^xx — Fyy — 2^2:y- 
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VI. ANALYSIS OF THE GROWTH EQUATIONS 



1, Ripple formation at high temperatures: 
Symmetric case 



This section is devoted to the study of the morpho- 
logical properties predicted by Eq. (40). This is not a 
simple task, due to large number of linear and nonlinear 
terms, each of which influence the surface morphology. 
The complexity of the problem is illustrated by some spe- 
cial cases of Eq. (^), for which the behavior is better 
understood. For example, when nonlinear terms and the 



noise are neglected {(^x = S,y 



Xy = 0,?/ - 0), 



Eq. ( |4C| ) reduces to a linear generalization of BH the- 
ory, which predicts ripple formation. It is also known 
that the isotropic KS equation, obtained by taking i^^ = 
Vy,Dxx = Dyy = Dxy/2, and = Xy, asymptotically 
predicts kinetic roughening, with morphology and expo- 
nents similar to those seen experimentally in ion sputter- 
ing p8|,p9|. For positive i^x and i^y, Eq. (Eol) reduces to 



the anisotropic KPZ equation, whose scalin g b ehavior is 
controlled by the sign of the product Xx ■ Xy |£3[ . Finally, 



recent integration by Rost and Krug |6£ 
anisotropic KS equation (i.e., when 77 = 



of the noiseless 
0) showed that 



when Xx ■ Xy < 0, ripples unaccounted for by the linear 
theory appear, their direction being rotated with respect 
to the ion direction. 

To predict the morphology of ion-sputtered surfaces, 
we need to gain a full understanding of the behavior pre- 
dicted by ( pO| ) in the physically relevant two dimensional 
case, going beyond the special cases. Help is provided 
by the recent numerical integration of (^0|) by Park et al. 
1 100 1 that indicates a clear separation in time of the linear 
and nonlinear behaviors. The results show that before a 
characteristic time tc has been reached, the morphology is 
fully described by the linear theory, as if nonlinear terms 
were not present. However, after tc the nonlinear terms 
completely determine the surface morphology. These re- 
sults offer a natural layout for our discussion. In section 
VIA we will limit our discussion to the linear theory. 
However, even in this case we have to distinguish four 
different cases, depending on whether the surface diffu- 
sion in the system is thermally generated or of the effec- 
tive type associated with th e ion ero sion process. Conse- 



quently, in Sections VI A 1 - VI A 2 , we discuss the high 



temperature case, when relaxation is by thermal surface 
diffusion, treating separately the symmetric {c — fi), and 
asymmetric (ct h) cases. Next we turn our attention 
to low temperature ripple formation, when surface relax- 
ation is dominated by erosion, and we again distinguish 
the sym metric and asym metric c a ses (Se cts. VI A 3 and 
VIA4D - Finally, Sections |VI B l| - |VIB4| are devoted to 



the effect of the nonlinear terms, addressing such impor- 
tant features as ripple stabilization, rotated ripples and 
kinetic roughening. 



A. Linear theory 



In this section we discuss the process of ripple forma- 
tion in the symmetric case a = ^, when the relaxation is 
by thermally activated surface diffusion. Thus we assume 
that the magnitude of the thermally activated surface dif- 
fusion coefficient, K , is much larger than Dxx, Dxy, Dyy, 
generated by the ion bombardment process. This is al- 
ways the case for high temperatures since K increases as 
{l/T)exp{-AE/kBT) with T, while ion induced effec- 
tive smoothing terms are independent of T. Dropping 
the nonlinear terms in Eq. (HQ), we obtain 



dh dh d'^h d'^h 



dt 



dx^ dxd'^y 



KV''h + Tj{x,y,t), (56) 



where the coefficients can be obtained from Eqs. ( pl| ) 
( ^ ) by taking a = fj,: 

vo = Fc, 

j = Fsiay-l), 

Fa 2 

= — , 

Fa'^s , 



^2 = ^s{aU3c'-3s^) + ays^~3}. 
bai 



(57) 



Since the surface morphology depends on the signs and 
absolute values of the coefficients in Eq. ( |56| ) , in the fol- 
lowing we discuss in detail their behavior as functions 
of the angle of incidence and the reduced penetration 
depth Ucr- 

Erosion velocity, vq: The vq term describes the ero- 
sion velocity of a flat surface. This term does not affect 
the ripple characteristics, such as ripple wavelength and 
ripple amplitude, and can be eliminated from the sur- 
face evolution equation by the coordinate transformation 
h — h + Vat. This corresponds to a transformation to the 
coordinate frame moving with the eroded surface. How- 
ever, since vo is the largest contribution to the erosion 
rate and is the only one that contributes in the linear 
theory, it is worthwhile to investigate its dependence on 
9 and a^. Fig. |^ shows the vq dependence on the angle 
of incidence 9 for three different values of the reduced 
penetration depth a^. From Eq. (|57|), vq is positive for 
all 9 and a^-. In experiments wo(^) corresponds to the 
secondary ion yield variation with the incidence angle 9, 
i.e. vq{9) = JYfiat{9)/n, where n is the density of target 
atoms. 
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FIG. 6. The erosion velocity, vo/H as a function of 6. 
The three curves correspond to the reduced penetration 
depth acr ~ 1 (solid line), Oo- = 2 (dotted line), ao- = 3 
(dashed line). The velocity has been normalized by a factor 
n — peJ/(\/27ra), independent of 9. 

Note that vq has the characteristic increasing part for 
small 9, followed by saturation and decrease for large 9, 
similar to the measured yield ||] . 

Traveling ripples, 7, Qi, Q2: If we consider a periodic 
perturbation with wave vector (q^r, qy) in the form 



h = —vqI + Aexp [i{qxX + qy-y — ujt) + rt] 



from Eq. (56) we obtain the mode velocity 
and the growth rate 



(58) 



(59) 



(60) 



Thus the coefficients 7,i7i,r22 contribute to the Fourier 
mode velocity uj in an anisotropic way that reflects the 
asymmetry of the x and y directions for oblique (9 ^ 0) 
ion incidence. The coefficients VxTVy^K, on the other 
hand, contribute to t he growth rate of the mode ampli- 
tude. Carter |^,101| pointed out that dispersive terms, 
such as rii and U,2, destroy the translational invariance 
of the periodic morphology because the different modes 
travel with different velocities. Note, however, that the 
existence of a ripple structure means that there is essen- 
tially only one Fourier mode describing the surface mor- 
phology, which will thus move across the surface with 
velocity to. The coefficient 7 contributes only to the ve- 
locity of the ripples along the x direction, leaving unaf- 
fected the y component of the ripple velocity. Thus, as 
expected, 7 = for normal incidence (6* = 0). Similarly 
to the Vq term, 7 does not affect the ripple characteristics 
and can actually be eliminated using the transformation 
h = h{x — 7t, t). 



FIG. 7. The coefficient 7/n as a function of the angle of in- 
cidence 9 for three reduced penetration depths: = 1 (solid 
line); = 2 (dotted line); ao- = 3 (dashed line). 

As can be seen in Fig. 0, 7 can change sign with 9, in- 
dicating that ripples travel in both positive and negative 
directions along the x coordinate, depending on the angle 
of incidence and the penetration depth: ripples travel in 
the positive x direction for small 6 and in the negative x 
direction for larger 9. Travelling ripples were observed in 
numerical simulations of Koponen et al. |]90| . 

As discussed above, the terms fii, Q.2 also contribute 
to the travelling of the ripples, and thus have no further 
effect on the surface morphology. Fig. ^ shows the coef- 
ficients Vti and 0,2 as functions of the angle of incidence 
9. We find that the absolute value of these coefficients at 
small angles is small compared to 7 (see Fig. |^), thus the 
main contribution to the ripple velocity comes from the 
["fdh/dx) term. On the other hand, for angles 9 > 60°, 
these terms are comparable to or larger than 7. 




FIG. 8. The reduced third order coefficients Q.\/Il (a) and 
Q.2 /n (b) as functions of the angle of incidence 9 for three re- 
duced penetration depths: ~ 2 (solid line); ~ S (dotted 
line); = 4 (dashed line). 
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FIG. 9. The surface tension coefficients Vx/Ii (a) and I'y/Ii 
(b) as functions of the angle of incidence 9 for three reduced 
penetration depths: Ocr = 2 (solid line); — 3 (dotted line); 
au — 4 (dashed line). 

The coefficients and Vy : As we discussed above (see 
Sect. IV B) the negative surface tension coefficients are 



the origin of the instabiUty responsible for ripple forma- 
tion. Consequently, they play a particularly important 
role in determining the surface morphology. The coeffi- 
cients Vx and Vy are not equal due to the fact that the 
direction of the ion beam breaks the symmetry along the 
surface. As seen in Eq. (57), Vy is always negative, while 
Vx can change sign as 6 and vary, as shown in Fig. |^. 
The sign and the magnitude of Ux and Vy determine both 
the wavelength and the orientation of the ripples. 

Ripple wavelength and orientation: The experimental 
studies on ripple formation have mainly focused on the 
measurement of the ripple characteristics, such as the rip- 
ple wavelength and amplitude. Thus, a successful theory 
must address and predict these quantities. In the follow- 
ing we outline the method for calculating the ripple wave- 
lengths ix and iy. Taking into account the noisy charac- 
ter of Eq. ( |56| ) , the experimentally observed ripple wave- 
length corresponds to the unstable Fourier mode which 
yields the maximum value of the structure factor. The 
structure factor, 5(q, i), is calculated from the Fourier 
transform h{(\, t) of the instantaneous surface profile and 
is defined as 



5(q,i) = (/i(-q,t)/^(q,i)). 



where 



/i(q, t) 



dr 



exp(iqr)/i(r, t). 



(61) 



(62) 



Fourier transforming Eq. ( p6| ) and inserting the expres- 
sion for the Fourier transforms of h{r,t) into (|6l|), we 
obtain 

^(q,i) = (M-q,t)Mq,^))--|^^^^^^^, (63) 

where r is the growth rate of the mode q given by Eq. 
( |60| ) and is positive for all unstable Fourier modes in the 



system. We find that, depending on the sign of Vx and 
the relative magnitude of Vx and Vy, we can distinguish 
two cases: 

(i) For Vx < '-'y < 0, which, according to Eq. ([5^), 
holds when 



(64) 



the ripple structure is oriented in the x direction, with 
ripple wavelength 



Ix = 27r 




(65) 



This means that the maximum of ^(q, t) is at {J 0). 



To illustrate this, in Fig. 10 we show the dependence of 
the structure factor on the wavevectors qx and qy. The 
contour plot indicates the existence of a global maxi- 
mum at {\J 0), indicating that the ripples are ori- 
ented along the x direction. 

(ii) For Vx > Vy, which holds when 



(66) 



the ripple structure is oriented along the y-direction, with 
ripple wavelength 



iy = 2TT 




(67) 



_>• 




FIG. 10. Contour plot of the structure factor 2S{qx,qy)/J 
as a function of the two dimensionless wavevectors q^a and 
qyU calculated for the angle of incidence 9 — 30°. The reduced 
coefficients Vx/^ and Uy/Tl are obtained using Eq. (^7[), their 
values being Vx/^ ~ —0.057, and Vy/Ii = —0.0418, while 
K/n is taken to be 0.01. These parameter values correspond 
to Region I in Fig. |l^. 
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FIG. 11. Contour plot of the structure factor 2S{qx,qy)/ J 
as a function of the two dimensionless wavevectors q^a and 
Qya calculated for the angle of incidence — 60°. The reduced 
coefficients Ux/Il and Uy/Tl are obtained using Eq. (p7[), their 
values being u^/H = 0.0758, and Uy/H = -0.0379, while K/H 
is taken to be 0.01. These parameter values correspond to Re- 
gion II in Fig. 

Figure |l^ shows an example of this regime, indicating 

the existence of a global maximum at point (0, ■^), 
corresponding to the ripple structure oriented along the 
y direction. 

Phase diagram for ripple orientation — The results ob- 
tained on ripple formation can be summarized in a (0, a^) 
morphological phase diagram, shown in Fig. |l2[ which 
has the following regions: 

Region I — For small both and Vy are negative 
such that i>x < i^y, thus the ripples are oriented along the 
X direction. Their wavelength is £x = 2-K\j2Kj\vx\. 
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FIG. 12. Ripple orientation phase diagram for the isotropic 
case a — — Region I: <Vy < 0; Region II: v^. > Vy 




region is defined by [o-a 



FIG. 13. Ripple wavelengths £x (solid line) and £y (dot- 
ted line) as functions of the angle of incidence 6 for 
KpeJ I {\plTi<i) = 0.01. The reduced penetration depth is 
taken as acr ~ 2. 

The amplitude of the ripples is expected to be weakly 
modulated by the larger wavelength £y = ^-n l\vy\. 
The ripple amplitude grows as ft-o exp(rj;i), where 

). The boundary of this 
Vy{aa, 6), i.e. 

(68) 



Region II — This region is characterized by Vy < Vx- 
The ripples are directed along the y direction and have 
wavelength £y = 2TT^j2K/\vy\. Note that Region II ex- 
tends down to small values of the incidence angle 9 for 
small enough reduced penetration depth a^. This some- 
what unphysical result is a consequence of the assump- 
tion of a symmetric (cr = /i) distribution of deposited 
energy. We will see in the next section that the more 
physical asymmetric case with a > fj, leads in most cases 
to ripples only oriented along the x direction for small 
enough angles of incidence, as generally observed. 

Figure ^ shows the 9 dependence of the ripple wave- 
lengths along the x and y directions. In the framework 
of this model, where thermal surface diffusion is the only 
smoothing mechanism, the observed ripple orientation 
corresponds to the direction featuring the smallest value 
of £, and changes when = £y The prediction for the 
ripple wavelength close to 90° is q uestionable since re- 
flection ^9 70 1 and shadowing |102|, not incorporated in 
the model, start to play an important role during ion- 
bombardment at these high angles. 

Summary: The dependence of £ on the main physical 
parameters characterizing the sputtering process is given 
by 
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(69) 



This prediction has a number of consequences, some 
of w hich have been verified experimentally (see Section 

(a) Si nce th e penetration depth, a, is proportional to 
e^™ (see |IV E] ), and ct ~ ~ a, we have ~ const, and 
F ~ {ea)/{afj.) ^ e^^^™, which is independent of e, when 
m = 1/2. Consequently 



-1/2 



(70) 



i.e. the ripple wavelength is expected to decrease with 
the ion energy. 

(b) Taking into account that K is independent of the 
flux and v ^ J, we obtain that the ripple wavelength is 
also a decreasing function of the incident ion flux, given 
by 



p8| , |69 70| T > 1 range. As one can observe the qualitative 
behavior of and Vy is similar to that observed in the 
symmetric case. One interesting feature, however, must 
be emphasized: the increasing asymmetry leads to larger 
ripple wavelength, since the absolute values of and Vy 
decrease. With respect to the the third order linear terms 
rii and fl2, their behavior as functions of the angle of in- 
cidence can be also seen to be qualitatively analogous 
to the symmetric case. Thus the asymmetry does not 
change our conclusions regarding the travelling ripples. 
Ripple wavelength: The calculation of ripple characteris- 
tics in the asymmetric case is identical to the one used in 
the symmetric case. Therefore, we limit ourselves to the 
presentation of the results. Again, there are two possi- 
ble ripple directions, and the dominant one can be found 
from the maximum of the structure factor ( |63| ) or, as 
can be seen to be equivalent, from the maximum of the 
growth rate ( |60| ) : 

(i) When < Vy < 0, i.e.. 



Jl/2- 



(71) 



(c) As was mentioned above, the negative surface ten- 
sion is the origin of the instability leading to ripple for- 
mation. When both and Vy are negative the experi- 
mentally observed ripple structure has the direction for 
which the growth rate r is largest. However, in general, 
we expect a superposition of both wavelengths, where 
the long wavelength will appear as a modulation of the 
ripple amplitude. Indeed, such modulations have been 



observed both experimentally and numerically |100|. 

(d) An important prediction of this model, illustrated 
in Fig. |lj, is the existence of the critical angle 6c where 
the ripple orientation changes. In the case when surface 
diffusion is thermally activated, this transition coincides 
with the condition v-r — Vn- 



/s2(2 + c2) +t2c2(1 + 2c2) -r4c4 , , 

aa> \l — , (72) 



the ripple structure is oriented along the a;-direction with 
ripple wavelength = 27ry^-^. 

(ii) When > i^y, i-e.. 



/s2(2 + c2) +t2c2(H-2c2) -t4c4 ^ ^ 

a<T < \l , (73) 



the ripples are oriented along the y-direction, with ripple 



wavelength £y = 27ri / 



2. Ripple formation at high temperatures: 
Asymmetric case 

The results of the previous section were derived for the 
isotropic case, cr = /x. While this approximation consider- 
ably simplifies our discussion, most systems present some 
anisotropy in the deposited energy distribution. In this 
section we demonstrate that the existence of anisotropy 
does not modify the overall qualitative result on the ex- 
istence of the two parameter regions corresponding to 
ripples oriented along the x or y directions. However, 
anisotropy does change the numerical value of the rip- 
ple wavelength and the exact boundary between the two 
morphological regions: we demonstrate that, for large 
enough anisotropics, if the incidence angle is small only 
ripples oriented along the x direction are possible. 

Fig. 14 shows the coefficients i>x and z/j,, given by Eq. 
(p7|), as functions of the angle of incidence 9, for three 
different degrees of asymmetry r = a/fi in the physical 
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FIG. 14. The effective surface tensions Vx/Tl (a) and z^y/II 
(b) plotted as functions of the incidence angle 6 in the asym- 
metric case. The curves correspond to values of the asymme- 
try parameter r — 1.5 (solid line), t — 2 (dotted line), t — 3 
(dashed line), with — 2. 
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Phase diagram for ripple orientation: To consider the 
effect of asymmetry on the different regimes in ion sput- 
tering, we have studied the ripple orientation phase dia- 
gram for different values of t. As t changes, we find a 
smooth evolution which docs not uncover any new mor- 
phological regime. However, the topology of the phase 
diagram does change as t increases. For r < \/3 the 
topology of the phase diagram is similar to the symmet- 
ric case (see Fig. ^2|). As Fig. |l^ illustrates, for r > \/3 
the ripples oriented along the y direction, predicted by 
the linear theory for small enough 9 and Cq-, are absent, 
which is consistent with most experimental observations. 
The characteristics of Region / and Region // of the 
phase diagram are the same as in the symmetric case. 
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FIG. 15. Ripple orientation phase diagram for the asym- 
metric case with asymmetry parameter r = 2 and = 2. 
Region I; Vx < Vy < 0; Region II; > Vy- 



believe this explains the rip ples observed at low temper- 
atures both experimentally |12(| and in computer simula- 
tions @. 

Neglecting the thermally induced relaxation terms (i.e, 
taking K — 0), nonlinear terms and the terms wq, 7, f^i 
and r22. that do not affect the ripple characteristics, from 
Eq. ( |40| ) we obtain the linear equation 



dh 
9t 



dx^ ^ dip 



-D 



xy 



(74) 



The expressions for the coefficients of the ion-induced ef- 
fective smoothing terms can be obtained from Eqs. (^l|)- 
(^) using a ~ fi: 



D 



yy 
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J {ats-^c^ + aliec^s^ - As^) + Sc^ - 123^} , 



6{a2s2c2+c2-2s2} 



(75) 



From Eq. (75), Dyy is always positive, while D^y and 



Dxx change sign with 9. As we discuss below, the posi- 
tive Dxx and Dyy coefficients play a role similar to ther- 
mally activated surface diffusion. Furthermore, the ab- 
solute value of Dxx is comparable with the thermally 
activated surface di ffusi on coefficient even at high tem- 
peratures (see Sect. VI] ). 

Ripple wavelength and orientation: The ripple wave- 
length and orientation can be calcul ated following the 
arguments presented in Section VI A 1 , being determined 
by the maxima of the structure factor S{q,t). The 
growth rate r is now given by 



3. Ripple formation at low temperatures: 
Symmetric case 

In the previous two sections we discussed the process of 
ripple formation when the origin of surface smoothing is 
surface diffusion, described by the —KW^h term. How- 
ever, in the series expansion of the erosion velocity we 
found linear fourth order terms of the form —D^xd'^h, 
—Dxyd^dyh, and —Dyydyh, which are formally equiv- 
alent to the thermally induced surface diffusion terms. 
These terms arise as a higher order correction to the lo- 
cal surface curvature, being fully determined by the pro- 
cess of surface erosion. Consequently, these terms do not 
imply actual mass transport along the surface, as ther- 
mal surface diffusion does. In this section we show that, 
in some parameter regions, these terms have a smooth- 
ing effect that counteracts the erosion instability, in such 
a way that they can also lead to ripple formation. We 



riqx,qy) = - {I'xqx + ^yQy + D^^q^ 



+Dxyqlql 



Dyyqt). 



(76) 



In principle the asymmetry of the Dij coefficients may 
lead to maxima of <S'(q, t) at nonzero qx and qy values, 
which correspond to ripples forming a nonzero angle with 
both the X and y directions. However, straightforward 
calculations indicate that the following condition holds 



(77) 



for all values of Go- and 9 in ( |57| ) and (|7^). This iden- 
tity implies that no extrema (q*, q*) of S'(q, t) exist other 
than of the form {q*, 0) or (0, q*). Of these two solutions 
the one with the largest positive value of r{qx, qy) corre- 
sponds to the observed ripple structure. For small angles 
of incidence so that Dxx > (Region I in Fig. p^, 
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FIG. 16. Morphological phase diagram in the symmetric 
case a = fi = 1. Different regions correspond to: Region I: 
i^x < 0, Uy < 0, > 0, Dyy > and > r^; Region II: 



ly-x <0,Uy< 0, < 0, Dyy > 0; Region III: 
Dxx < and Dyy > 0. 



>0, 



<0, 



it can be easily seen that < 0, and the absolute max- 
imum of 5(q, t) is at (g* , 0) with q* = y^li'xl/'^Dxx, thus 
the ripple structure is aligned along the x direction (see 
for example Fig. p7[ ). 




FIG. 18. Contour plot of the structure factor 2S{qx,qy)/J 
shown as a function of the two dimensionless wavevectors q^a 
and qya calculated for the angle of incidence 6 = 30°. The re- 
duced coefficients Vx/^ and Vy/U are obtained using Eq. (p7|), 
their values being u^/Ii = —0.0758 and /^y/II = —0.0379. 
These parameters values correspond to Region II in Fig. 



Crossing the D^x = line into Region II in Fig. 16 the 
(g*, 0) solution disappears, the structure factor having an 
extremum at {0,q*), with q* = ^J\vy\l2Dyy. However, 
this extremum is not a global maximum, see for example 
Fig. |l^. The lower boundary of Region II is provided by 
the condition i^x = 0. When we cross it (entering Region 
III in Fig. we have Dxx < and h'x > 0. Under this 
condition, again, there exists an extremum of ^(q, t) at 
(g*,0), with q* = ^ Vxl'2\Dxx\- However, the structure 
factor takes its absolute maximum at (0, g*). 

Phase diagram for ripple orientation: In summary, 
three different regions can be determined in the morpho- 
logical phase diagram shown in Fig. |l^ for the case of 
ion-induced effective smoothing in the symmetric a — fi 
case. 

Region I — In this region, the surface tension coeffi- 
cients Vx and Vy are negative, while Dxx and Dyy are 
both positive. The observed ripple structure corresponds 
to the maximum of S'(q, t), which indicates that the rip- 
ple structure is oriented along the x direction. The lower 
boundary of this region, separating it from Region II, is 
given by the Dxx^cLa^O) = line or, equivalently, by 



(2s2 - 3c2) + V6c4 + 4s4 



(78) 



FIG. 17. Contour plot of the structure factor 2S{qx,qy)/ J 
shown as a function of the two dimensionless wavevectors qxa 
and Qya calculated for the angle of incidence 6 = 30°. The re- 
duced coefficients Vx/Tl and I'y/Tl are obtained using Eq. (^7|), 
their values being Ux/U = -0.0578 and Uy/R = -0.0418. 
These parameter values correspond to Region I in Fig. |l^. 



Region II — In this region both Dxx and Vx are nega- 
tive. This region is bounded below by the Vx{a„^6) = 
line, given by 



(2s2 - c2) 



(79) 



In a continuum description, the maximum of r(q) is at 
infinite q, thus our theory possibly breaks down in the 
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sense that not even non-linear effects can be expected 
to stabilize the surface under such conditions. In such a 
case, a higher order Taylor expansion should be carried 
out in Sec. ^ in order to be able to describe our system. 
Additional effects, such as shadowing, could also start to 
play a role under such conditions. 

Region III — In this region D^x is negative and is 
positive. Thus the instability given by the negative D^x 
is smoothed out by the positive Vx- Since the structure 
factor takes on its absolute maximum at the finite wave 
vector (0, q*), in principle there is still a ripple structure 
oriented along the y direction. However, remarks similar 
to those made in Region II might apply here, since we 
still have Dxx < 0. 

Summary: In the presence of ion-induced effective 
smoothing the dependence of the ripple wavelength on 
the physical parameters is diffe rent fro m the case of ther- 
mal surface diffusion (see Sect. VI A 1 ) . Here we summa- 
rize some of the differences. 

(a) The dependence of i on the ion energy is given by 



(80) 



indicating that the ripple wavelength depends linearly on 
the penetration depth a. This is very different from the 
behavior predicted by Eq. (|5^), derived for thermal sur- 
face diffusion, and indicates that monitoring the ripple 
wavelength dependence on e can be used to identify the 
dominant smoothing mechanism. Such a linear behavior 
of £ on e has indeed been seen experimentally (see Section 




[HAp . 

(b) From Eq. ( ^ ) it also follows that the ripple wave- 
length is independent of the incident ion flux. This pre- 
diction is again quite different from the case dominated 
by thermal surface diffusion, given by Eq. (|7l|). Such a 
flux independent behavior has been observed experimen- 
tally (see Section [HA). 



Phase diagram for ripple orientation — The topology 
of the morphological phase diagram and the character- 
istics of the three main regions are not changed by the 
asymmetry. We find that the only effect of the asymme- 
try is to move the boundaries smoothly to larger values 
of 6* as r increases. The condition D{aa-, 9) — (see Eq. 
([ZS])) now takes the form 



+ ^ ^2) -^3-5 , 



(81) 



and the condition Vx{a„,9) = (see Eq. (79)) becomes 

(82) 



2s4 ^2^2g2 _ ^4^4 



Finally, anal ogues o f characteristics (c) and (d), dis- 
cussed in Sect. VI A 1, apply here as well. 



B. Nonlinear theory 

As we demonstrated in the previous section, linear the- 
ory can predict many features of ripple formation, such 
as the ripple wavelength and orientation, both at high 
and low temperatures. However, a number of impor- 
tant experimental features are incorrectly predicted by 
linear theory. They include the stabilization of the rip- 
ple amplitude (according to the linear theory the ampli- 
tude increases indefinitely at an exponential rate) or the 
presence of kinetic roughening (completely unaccounted 
for by the linear approach). To account for these fea- 
tures, we have to consider the effect of the nonlinear 
terms. Consequently, this section is devoted to the ef- 
fect of the nonlinear terms on the surface morphology. 
There is an important difference between the linear and 
nonlinear theories: while all predictions of the linear the- 
ory can be calculated analytically (as we demonstrated 
in the previous section), the discussion of the nonlinear 
effects requires a combination of analytical and numer- 
ical tools. Even with these tools, the understanding of 
the nonlinear effects is far less complete than that of the 
linear theory. 



4. Ripple formation at low temperatures: 
Asymmetric case 

In this section we discuss the effect of asymmetry 
((T 7^ /i) of the energy distribution on the morphological 
regimes predicted by Eq. ([74|). Wc find that the coeffi- 
cients of Eq. (jrj) vary slowly with the asymmetry, but 
this does not change the qualitative picture presented in 
the previous section, regarding the ripple wavelength and 
orientation, or the major morphological regimes found in 
the isotropic case, including the conditions of validity 
of our continuum theory. Specifically, we find that the 
asymmetry enlarges the region in 9 where Dxx and Dyy 
are positive, thus shifting Region // to larger values of 



1. High temperature morphology: Symmetric case 

In the high temperature regime, where thermal surface 
diffusion dominates over ion-induced effective smoothing, 
the nonlinear equation of the interface evolution has the 
form 



dh dh ^ 

(dh\ fd'^h 



9x / \dy 
d^h Xx (dh 



dh\ fd^h 
dx ) \ dx'^ 
d^h 



\y ( dh 



dy 



+ni^ + n2^^^~KV^h + Tj{x,y,t), (83) 
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where the coefRcients of the hnear terms, vq, 7, tz-m 
rii, r^2, and K have been discussed in Sections VI A 1 



and VI A 2, The coefRcients of the nonlinear terms in 



the symmetric case (ct = fi) are 
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(84) 
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Next we discuss the physical interpretation and the 
behavior of these coefficients as functions of 9 and a^. 

The coefficients o,nd ^y: Fig. |l^ shows the nonlin- 
ear coefficients (^x and as functions of the angle of 
incidence 9. As the numerical analysis of Eq. (|8^) shows, 
these nonlinearities are responsible for the development 
of overhangs on the surface |103[ . Even though the S^x 
and terms are expected not to determine the asymp- 
totic scaling behavior, they can play a relevant role at 
intermediate time scales after the development of the rip- 
ple structure, particularly in the region of large 9. The 
precise contribution of these nonlinearities to t he surface 
dynamics is currently under investigation |l03t| - 

The coefficients Xx and Aj,; As we discussed in Sec. 



IV A 2, the morphology of the surface described by Eq. 



(^3|) depends on the relative signs of the the nonlinear 
terms and Aj^. As it is evident from Eq. (^4|), Aj, 
is negative for all angles of incidence and penetration 
depths. 
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FIG. 19. The reduced coefficients ^^/H (a) and ^y/W (b) 
shown as functions of 6 for different values of the reduced 
penetration depths: = 2 (solid line); = 3 (dotted line); 
aa = 4 (dashed line). 



FIG. 20. The reduced coefficients K/Tl (a) and A„/n (b) 
shown as functions of 9. The different curves correspond to 
different values of the reduced penetration depth: aa = 2 
(solid line); flo- = 3 (dotted line); ao- = 4 (dashed line). 

However, as shown in Fig. |2^, the sign of A^; depends 
on 9 and Oo-: A^; is negative for small 9 and changes sign 
for larger angles of incidence. In principle the nonlin- 
ear terms completely determine the surface morphology. 
Since the nonlinear terms are always present, an impor- 
tant question is whether the linear r egim es are relevant 
at all. Recent results by Park et al. [100| indicate that, 
while the nonlinear effects indeed change the surface mor- 
phology, the regime described by the linear terms is still 
visible for a wide range of parameters. By numerically in- 
tegrating Eq. ( p3| ) they have shown that there is a clear 
separation of the linear and nonlinear regimes in time: 
for times up to a crossover time tc the surface erodes 
as if the nonlinear terms would be completely absent, 
following the predictions of the linear theory. After tc, 
however, the nonlinear terms take over and completely 
determine the surface morphology. The transition from 
the linear to the nonlinear regime can be seen either by 
monitoring the surface width (which is proportional to 
the ripple amplitude) or the erosion velocity. The sim- 
ulations indicate that the width increases exponentially 
with time, as predicted by the linear theory, until t^, af- 



ter which the width growths at a much slower rate [100|. 
This transition is typically accompanied by the disap- 
pearance of ripples predicted by the linear theory and 
the appearance of either kinetic roughening or of a new 
rotated ripple structure. The erosion velocity is constant 
in the linear regime (before tc), while it increases or de- 
creases after t^ depending on the relative signs of the 
nonlinear terms. 



Crossover time: The crossover time t^. from the lin- 
ear to the nonlinear behavior can be estimated |10C] by 
comparing the strength of the linear terms with that of 
the nonlinear terms. Let the typical surface width at the 
crossover time tc be Wo — \/W'^{L, tc). Then from the 
linear equation we obtain 
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FIG. 21. Phase diagram for the isotropic case a — fj, — 1. 
Region I: < 0, fy < 0, Xx < 0, \y < 0; Region II: i/^ < 0, 
Uy < 0, > 0, Xy < 0; Region III: > 0, Vy < 0, X^ > 0, 

Xy < 0. 

Wo ~ eM'^tc/f ), (85) 
while from dth ~ A(V/i)^ we estimate 

Wo/t, ^ XWlie. (86) 



Combining these relations we obtain 



(87) 



In this expression, ly, K, and A correspond to the direc- 
tion parallel to the ripple orientation. The predicted A 
dependence of tc has been confirmed by numerical simu- 



lations 1 100 1 



Surface morphology: The surface morphology in the 
nonlinear regime depends on the relative signs of v^, I'y, 
Xx and Aj,. The different morphological regimes can be 
summarized in a phase diagram, shown in Fig. Next 
we discuss each of the phases separately. 

Region I — For small the nonlinear terms A^; and Aj, 
have the same (negative) sign, the boundary of this re- 
gion being given by the condition that Xx{aa,0) — 0, or 
equivalently 



(88) 



In this region both and Vy are negative, thus at shor t 
time scales {t < tc), the linear theory (see Section VI A 1 ) 
predicts ripples oriented along the direction (x for large 
Go- and y for small Ua-) for which the absolute value of v 

is largest, with ripple wavelength £ = 2tt^^. On the 

other hand, at long times {t 3> tc), the ripple structure 
disappears and the surface undergoes kinetic roughen- 



ing [100|. Since Aa; • Ay > 0, we expect the dynamics 



KPZ equation, i.e. the surface width follows W ^ L", 
W ^ t^, where the scaling exponents are a ~ 0.38 and 
(3 ~ 0.25 (see Sect. [IVAlD . 

Region II — In this region both the and the Vy co- 
efficients are still negative, but in contrast with Region I 
the product X^ ■ Xy is negative. Recent studies by Park 



et al. I IOC] have shown that in time three morphological 
regimes can be distinguished. For short times, t < tj, the 
rippl e struc ture predicted by the linear theory (see Sec- 
tion [V A 1 ) is observed, with ripples oriented along the 

For intcrme- 
If this 



of the kinetic roughening regime to be described by the 



direction which has the largest value of 
diate times tj < t < tjj, the surface is rough 
roughness would follow kinetic roughening, one would 
expect logarithmic scaling, described by the Edwards- 
Wilkinson equation (see Sections VI B 2- VI B 4), since 
Aa; ■ Ay < 0. However, this transient regime is somewhat 
different from what we expect during kinetic roughening. 
The numerical simulations often show the development 
of individual ripples, which soon disappear, and no long- 
range order is present in the system. However, at a sec- 
ond crossover time, i//, a new ripple structure suddenly 
forms, in which the ripples are stable and rotated an an- 

t6c with respect to the x direction. Rost and Krug 
have shown [for the deterministic limit of Eq. (|8^)] 
that, by defining = v^jvy and /3a = Xx/Xy, the fact 
that Qfi, > and fix < Q throughout Region II implies 
the existence of "cancellation modes" which dominate 
the dynamics and lead to this rotated ripple morphol- 
ogy. (Note the parameter ratios an and (3\ are not to 
be confused with the roughness and growth exponents a 
and [3 introduced in Section ||.) The angle 9c calculated 
by Rost and Krug has the value 6c — tan~^ \/—Xx/Xy 
(see also Appendix ^ , and can be obtained by moving 
to a rotated frame of coordinates that cancels the non- 
linear terms in the transverse direction. The boundary 
of Region // is given by the condition i^xio-a, 0) = 0, Eq. 

Region III — This region is characterized by a positive 
Ux and a negative Vy. At short time scales, t < tc, the pe- 
riodic structure associated with the instability is oriented 

along the y direction and has wavelength £y = 

Much less is known, however, about the nonlinear regime. 
Such an anisotropic and linearly unstable equation is un- 
explored in the context of growth equations. The analysis 
by Rost and Krug for the corresponding determinis- 
tic equation predicts that, given that (3x < Ui, < Q does 
hold all over Region HI, again cancellation modes induce 
a rotated ripple morphology. 

Summary: Even though several aspects of the scaling 
behavior predicted by Eq. ( ^3| ) remain to be clarified, we 
believe that this equation contains the relevant ingredi- 
ents for understanding roughening by ion bombardment. 
To summarize, at short time scales the morphology con- 
sists of a periodic structure oriented along the direction 
determined by the largest in absolute value of the nega- 
tive surface tension coefficients 
of Co- or changes the orientation of the ripples 



. Modifying the values 
At 
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long time scales we expect two different morphological 
regimes. One is characterized by the KPZ exponents, 
which might be observed in Region / in Fig. ^ In- 
deed, the values of the exponents reported by Eklund et 
al. |2^,|2^ are consistent within the experimental errors 
with the KPZ exponents. In Region // kinetic roughen- 
ing is not expected. Rather, nonlinear terms lead to a 
new ripple structure that is rotated with respect to the 
ion direction. Region III is less understood; analysis of 
the deterministic equation [6^ ] again predicts a rotated 
ripple structure. By tuning the values of 9 and/or a^, 
one may induce transitions among these morphological 
regimes. 



2. High temperature morphology: Asymmetric case 

In this section we discuss the effect of asymmetry on 
the scaling regimes of Eq. (^3|). Here again we obtain 
that the effect of asymmetry does not bring in new qual- 
itative features. Specifically, we find that the qualitative 
behavior of S^x and S^y is not affected by the asymmetry. 
As the asymmetry grows, the absolute value of the co- 
efficients in the region of small angles decreases and the 
peaks at large 9 increase. Similarly, while the numerical 
values of A^; and Ay are affected by the asymmetry t, 
their qualitative features are not. 

Finally, we find that the morphological diagram is 
topologically equivalent to the phase diagram obtained 
in the symmetric case (see Fig. |2^), the only difference 
being in the position of the boundaries. As r changes, we 
find a smooth evolution of the boundaries, which does not 
uncover any new regimes. Since the morphological prop- 
erties of the system in the three regimes are the same as 
those discussed in the symmetric case, we will not discuss 
them again. 



3. Low temperature morphology: Symmetric case 

In this section we discuss the effect of the effective sur- 
face smoothing on the surface morphology in the nonlin- 
ear regime. In the absence of thermally activated surface 
diffusion, from Eq. (40) we obtain the following equation 
governing the morphology evolution 
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FIG. 22. Nonlinear phase diagram for the isotropic case 
a — jj, — 1 at low temperatures. Region L < 0, < 0, 
Dxx > 0, Dyy > and < 0, < 0; Region Ila: < 0, 
Uy < 0, Dxx < 0, Dyy > and A^; < 0, Xy < 0; Region lib: 
< 0, i^y < 0, Dxx < 0, Dyy > and A^; > 0, Ay < 0; Region 

III: l^x > 0, < 0, DxX < 0, Dyy > lUid Xx > 0, Xy < 0. 



The terms 7, u^, Vy, fii, ^2, ^x, ^y, K 



Ay, as well as 



xy 



the ion- induced effective smoothing coefficients D^xi D 
and Dyy have been discussed in the previous sections. 

In the following we discuss the morphological phase di- 
agram predicted by Eq. (^ and shown in Fig. |2^. The 
different regimes have the following characteristics: 

Region I: The surface tensions Vx and h'y are negative 
while Dxx and Dyy are positive, and A^;, Ay are both 
negativ e. This regime has been previously described in 
Section VI B 1 (Regime I in Fig. ^), the only difference 
being that here the ion-induced effective surface smooth- 
ing plays the role of K . The boundary of this region is 
given by i^,,(a,, 9) = 0, Eq. (^. 

Region Ila: Here Vx, i^yi Xx, and Ay are still negative, 
Dyy is positive, while Dxx < 0. Consequently, along the 
X direction the surface is unstable, all modes growing 
exponentially. However, the nonlinear terms A^; and Ay 
have the same sign. The nonlinear regime in this param- 
eter region has not yet been explored numerically, thus 
its morphology is unknown. The boundary of this region 
is given by Xx{acr,9) = 0, Eq. (H). 

is negative, Dyy is pos- 



Region lib: In this region Dx 



itive, Vx < 0, f„ < 0, and A^ > 0, Ay < 0. Thus, the only 



difference of this region with respect to Region Ha is that 
the nonlinear terms have different signs. Similarly to Re- 
gion Ila, nothing is known about the nonlinear behavior. 
The boundary of this region is given by i>x{o,aTS) = 0, 
Eq. (^). 

Region III: Here we have Vx > 0, Uy < 0, Dxx < 0, 
Dyy > 0, A2: > and Ay < 0. This region has similar 
features to Region III in the phase diagram of Fig. |2^, 
except for the negative value of the Dxx coefficient. 
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FIG. 23. Nonlinear phase diagram for the anisotropic case 
with r = 2 and = 2. Different regions in the diagram 
correspond to: Region I: < 0, i^y < 0, D^x > 0, Dyy > 
and Xx < 0, \y < 0; Region Ila: Ux < 0, Uy < 0, Dxx < 0, 
Dyy > and A^; < 0, \y < 0; Region lib: ^■j; < 0, i^y < 0, 
< 0, Dyy > and Aa; > 0, Ay < Oj Region lie: z^a; < 0, 
Uy < 0, Dxx > 0, Dyy > and A^, > 0, Ay < 0; Region III: 
i^x > 0, uy < 0, Dxx < 0, Dyy > and A^, > 0, Ay < 0. 



4- Low temperature morphology: Asymmetric case 

In this section we discuss the effect of asymmetry on 
the long distance properties of Eq. ( |89| ) . The effect of the 
asymmetry on the coefficients appearing in the equation 
were discussed earlier. Therefore we concentrate here 
merely on the morphological phase diagram predicted 
by Eq. (p9[ ) for the asymmetric case, which is displayed 
in Fig. As before, asymmetry (r ^ 1) leads to a 
smooth shift of the boundaries of the regions provided 
by the lines where the coefficients OxxiaajO), Vx{ac,9), 
and \x {o-a , 9) change sign. In the presence of effective 
smoothing, however, asymmetry in the deposited energy 
distribution induces the appearance of a fifth morpho- 
logical regime. This is caused by the smooth motion of 
the boundary determined by Xx{ac!,9) = 0, which in- 
tersects for some value of r the boundary defined by 
Dxx{a<j,9) = 0. Comparison of Fig. |23| and Fig. ^ il- 
lustrates how the boundaries move with r. Regions I, 
Ila, lib and III are analogous of those shown in Fig. |2^. 
Region lie in the phase diagram, on the other hand, is 
analogous of Region II of the high temperature phase di- 
agram, shown in Fig. ^l], and all the conclusions obtained 
in that section regarding the morphological properties in 
this regime also apply here. 



VII. COMPARISON WITH EXPERIMENTS 

In this section we compare the predictions of the the- 
ory presented in this paper with experimental results 



on ripple formation and surface roughening. For a bet- 
ter presentation, we choose to structure the material 
around well known features of the morphological evolu- 
tion, present the theoretical predictions and discuss to 
which extent are they supported by the available exper- 
imental data. We also discuss predictions that have not 
been tested in sufficient detail but could offer future tests 
of the theory. 

Ripple amplitude: A key quantity in ripple formation 
is the time evolution o f the r ipple amplitude. As we 
have shown in Section VIAl, at early times {t < tc) 
the ripple amplitude grows exponentially, following h ~ 
exp (r(g* , where r is the growth rate of the most 

unstable mode (qxily)- According to the linear theory, 
this growth continues indefinitely. In contrast, the non- 
linear theory predicts that the amplitude should stabilize 
after time where tc is given by Eq. (|87|). This is consis- 
tent with the experimental investigations [|l8|,^ . On the 
other hand, recently Erlebacher et al. also found that 
at initial stages the ripple morphology is growing expo- 
nentially. Furthermore, they observed that at some time 
tc the exponential growth stops and the ripple amplitude 
saturates. Measuring the temperature dependence of the 
saturation curves, they found that rescaling the time t 
with a factor i/^/AK and the amplitude h with yV/2i?, 
the different curves representing the amplitude as a func- 
tion of time collapse onto a single one. This result is in 
excellent agreement with our prediction that suggest that 
plotting the result in terms of the rescaled parameters, 
t/tc and h^y v/2K, the different curves should collapse 
100 1 . They also offer direct proof that the nonlinear 



terms play a major role in determining the amplitude of 
the ripples, indicating that the incorporation of the non- 
linear mechanisms in the theory of ripple formation is 
essential. 

Temperature dependence of the ripple wavelength: A 
key quantity that provides direct information about the 
nature of the relaxation mechanism is the temperature 
dependence of the ripple wavelength. Our results in- 
dicate that there are two mechanisms contributing to 
surface relaxati on: t hermally induced surface diffusion 



(Sects. VI A 1 - VI A 2 ) and ion- induced smoothing (Sects. 
VI A 3 - VI A 4 ) . At high temperatures thermal surface dif- 
fusion is rather intensive, thus it is the main mechanism 
determining the relaxation process, the ripple wavelength 
being given by (|6^) or (^. Since the surface diffusion 
constant K decreases with T as (1/T) exp(— AiJ/fc^r), 
the ripple wavelength is also expected to decrease expo- 
nentially with T . Indeed, such an exponential tempera- 
ture dependence of £ has bee n obser ved by various groups 
|]l2| , ^ . However, in Section VI A3 we demonstrated the 
existence of ion induced smoothing, that is present at 
any temperature. Thus, up to some inessential numer- 
ical factors, the total surface diffusion constants have a 
form = K + D, where D is independent of temper- 
ature. Since K decreases exponentially with T, at low 
enough temperatures we have K ^ D, indicating that 
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the main relaxation mechanism is ion-induced. Conse- 
quently, below a certain critical temperature Tc, given 
by K{Tc) — D, one expects the ripple wavelength to be 
independent of T. Support for this scenario has been 
provided by the experiments in [|l2| and and the 
molecular dynamics simulations of Koponen et al. [po| . 
Consequently, as the theoretical results in this paper in- 
dicate, the temperature dependence of i can be used to 
identify the relaxation mechanism: when t increases ex- 
ponentially with T, we are dealing with relaxation by 
thermal surface diffusion, while a temperature indepen- 
dent wavelength is an indication of ion induced smooth- 
ing. 

Ripple orientation: An important feature of ripple for- 
mation is that, as the linear theory predicts, the ripple 
orientation depends on the angle of incidence 9. The de- 
pendence of the ripple orientation on the experimental 
parameters has been summarized in the phase diagrams 
shown in Figs. ||^, H and ||. In general, 

for physical values of the asymmetry parameter {i.e. for 
T > 1), we find that for small angles the ripples are ori- 
ented in the x direction (along the incoming ions), and 
they change orientation to the y direction for large 9. The 
boundary separating these two morphological regions de- 
pends on the parameters characterizing the sputtering 
process, such as the ion penetration depth and the geome- 
try of the deposited energy distribution. Such transition 
in the ripple orientation has been found in the simula- 
tions of JSSl l, where for 9 < 45° the observed ripples were 
oriented along the x direction, while for 9 > 45°, they 
changed their orientation to y. Furthermore, the nonlin- 
ear theory predicts that after the nonlinear terms take 
over, new ripples, with orientation differen t from both 
X and y directions, might appear (see Sect. VI B 1 ) . To 
the best of our knowledge, such rotated ripples have not 
been observed experimentally as yet. Nevertheless, this 
morphology might also lead to additional effects, such as 
shadowing, which have been neglected in our approach. 

Ripple wavelength dependence on the flux: Depending 
on the nature of the relaxation mechanism, the linear 
theory has two different predictions on the flux depen- 
dence of the ripple wavelength: for high temperatures, 
when thermal surface diffusion dominates, one expects 
£ ^ J~^/^ [see Eq. (|7l])], while at low temperatures, char- 
acterized by ion induced smoothing, we expect the wave- 
length to be independent of flux (see Sect. VI A 3| ) . Con- 
sequently, due to its strong dependence on the relaxation 
mechanism, the flux dependence of the ripple wavelength 
can also be used to identify the relaxation mechanism. 
Indeed, a number of experiments (l^Jl^ are compatible 
with the prediction of a flux independent wavelength. 
Other aspects of ripple characteristics (such as energy or 
temperature dependence) also lead to the same conclu- 
sion. On the other hand, we are not aware of results 
indicating decreasing ripple wavelength with increasing 
flux. However, support for the relevance of thermal sur- 
face diffusion comes from the experiments of Chason et 
al. [^, who reported that the growth rate r{q*,q*) as 



a function of flux follows the predictions of the linear 
theory with thermal surface diffusion. 

Ripple wavelength dependence on the ion energy: The 
linear theory indicates that the ion energy dependence of 
the ripple wavelength can be used to distinguish between 
the two relaxation mechanisms: at high temperatures 
we expect £ ~ e~^/^ [see Eq. ([70|)], i.e., the wavelength 
decreases with the energy, while at low temperature we 
have £ e^™ ~ e [see Eq. (|80|)], i.e., the wavelength 
should increase with energy, strikingly different predic- 
tions. A number of experimental groups have found that 
the ripple wavelength increases linearly with the ion en- 
ergy [|7|,p|Jl^]. However, while we are not aware of any 
direct observation of a decreasing ripple wavelength with 
increasing ion energy, the growth rate dependence on the 
ion energy measured by Chason et al. pj^ provided re- 
sults which are in agreement with the predictions based 
on thermal surface diffusion. 

The magnitude of the effective surface diffusion con- 
stant: Since the transition between the low and high tem- 
perature regimes is determined by the relative magnitude 
oi D^x, D^y, Dyy {lon Induccd smoothing), and K (ther- 
mal surface diffusion), we need to estimate the magnitude 
of these constants. In the following we give an order of 
magnitude estimate for the effective surface diffusion con- 
stant Dxx and compare it to K, using data from [|l8|jl9[ 
for Si(OOl). Taking Y = 2.6, J = 670 ^A/cm^, e = 9 
keV, a = 100 A, a„ = 2, = 4, and 9 = 40°, Eq. 
( pl| ) gives Dxx — 12 X 10~^^ cm^/s. For comparison, 
at T = 550 C it is estimated ^ that ~ 34 x lO'^s 
cm^/s. Hence, since K decreases exponentially with tem- 
perature, ion induced smoothing can be significant at 
low temperatures (including room temperature) , in some 
cases being comparable or more relevant than thermal 
surface diffusion. 

Kinetic roughening: An important feature of our the- 
ory is that it goes beyond the linear approach, han- 
dling systematically th e non linear effects as well. As we 
demonstrated in Sect. VI B , the presence of the nonlin- 
ear terms can affect both the dynamics and the mor- 
phology of the surface. The first and the most dramatic 
consequence is the stabilization of the ripple amplitude, 
discussed above. Furthermore, after the stabilization of 
the ripple amplitude the surface morphology is rather 
different from the morphology predicted by the linear 
theory. In particular, depending on the signs of A^; and 
Aj,, different morphological features can develop. When 
Xx^y is positive, at large times the surface undergoes ki- 
netic roughening, following the predictions of the KPZ 
equation. This behavior has, indeed, been observed ex- 
perimentally 1 28 2^ and numerically [ ^ , providing di- 
rect support for the predictions of the nonlinear theory. 
When XxXy is negati ve, direct numerical integration of 



the nonlinear theory [ IOC ,|65| indicated the existence of a 
new, rotated ripple structure. The absence of experimen- 
tal data on this phase might be due to the required large 
sputtering times: the simulations indicate [100| that be- 
tween the linear regime and the formation of the rotated 
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ripple structure there is a rather long transient regime 
with an apparently rough surface morphology. The above 
predictions apply when the surface diffusion terms, ion 
or thermally induced, act to smooth the surface (i.e., 
D"^ > 0). However, at low temperatures, when ion in- 
duced smoothing dominates, surface diffusion can gener- 
ate an ins tability that can further modify this behavior 
(see Sect. VI B 3| ) . In general, while rather detailed ex- 
perimental data are available describing the linear regime 
of ion sputtering, explanation of the nonlinear regime is 
only at its beginning, hiding the possibility of new inter- 
esting phases and behaviors. 



VIII. CONCLUSIONS 



In this paper we investigated the morphological prop- 
erties of surfaces eroded by ion bombardment. Start- 
ing from the expression for the erosion velocity derived 
in the framework of Sigmund's theory of sputtering of 
amorphous targets, we derived a stochastic partial dif- 
ferential equation for the surface height, which involves 
up to fourth order derivatives of the height, and incor- 
porates surface diffusion and the fluctuations arising in 
the erosion process due to the inhomogeneities in the ion 
flux. In some special cases the derived nonlinear theory 
reduces to the much studied KS or the KPZ equations, 
well known descriptions of dynamically evolving surfaces. 
However, in contrast with these theories, which have been 
derived using symmetry and conservation considerations 
p7| , here we derived the continuum theory directly from 
a microscopic model of sputtering, and thus all coeffi- 
cients can be explicitly expressed in terms of the physi- 
cal parameters (such as angle of incidence, ion penetra- 
tion depth, etc.) characterizing the ion bombardment 
process. An important feature of the derived nonlinear 
continuum theory is that the linear and the nonlinear 
regimes are separated in time. As a result, they can be 
discussed separately, the former controlling the behavior 
at early times, the latter at late times. Furthermore, an 
important result of our calculations is that higher order 
effects of the sputtering process can smooth the surface. 
This effective mechanism was necessary to explain ripple 
formation at low temperatures, when thermally induced 
surface diffusion is not relevant. Consequently, based on 
these two ingredients (separation of time scales between 
linear and the nonlinear regimes and the existence of two 
different relaxation mechanisms) we have discussed four 
different cases. In the linear high temperature regime 
the equations reduce to the linear theory of Bradley and 
Harper, predicting ripple formation, and explaining such 
experimentally observed phenomenon as ripple orienta- 
tion (and its change with and 9) , exponential increase 
in ripple amplitude (valid for short times), or flux and en- 
ergy dependence of the ripple wavelength. On the other 
hand, phenomena not explained by this approach, such 
as the stabilization of the ripple amplitude, can be ex- 



plained by considering the nonlinear terms as well. We 
also show that, depending on the sign of the coefficients 
of the nonlinear terms, the late time morphology of the 
surface is either rough, or dominated by new rotated rip- 
ples. The rough phase is expected to be described by 
the KPZ equation, which has its own significance: while 
the introduction of the KPZ equation has catalyzed an 
explosion in the study of the morphological properties 
of growing surfaces, there are very few actual surfaces 
that are described by it (and not by one of its offsprings, 
such as the MBE or related equations [^), and most no- 
tably none, as far as we know, that describe crystal sur- 
faces. The sputtering problem provides one of the first 
systems that is convincingly described by this continuum 
theory. Many of the previous mysteries of low tempera- 
ture ripple formation have also been solved by the present 
theory. By deriving the higher order ion-induced effec- 
tive smoothing terms, we can explain the existence of 
ripples at temperatures where thermally induced surface 
diffusion is not active. We showed that the derived ef- 
fective smoothing affects both the linear and the nonlin- 
ear regimes, governing the early time ripple formation, 
and the late time nonlinear behavior. The coexistence of 
thermal and ion induced smoothing can explain the sta- 
bilization of the ripple wavelength at low temperatures, 
in contrast with its exponential T dependence at higher 
temperatures. On the other hand, our theory has limi- 
tations, most of which can be already identified in Sig- 
mund's and Bradley and Harper's theories, with which 
it is related. Namely, it is devised for amorphous sub- 
strates, whereupon it neglects effects such as viscous re- 
laxation 1?^] , which might be the cause for the failure of 
the theory to predict the absence in many experiments 
of ripples at low (but non-zero) angles of incidence. This 
issue should constitute one of the most important exten- 
sions of our present theory. Perhaps related with this, we 
have seen that there exist parameter regions at low tem- 
peratures within which our theory breaks down, due to 
the unstable higher order derivative terms that occur. A 
relevant issue is thus to determine the correct continuum 
description of the surface under these conditions. 



Most of the predictions offered by the presented con- 
tinuum theory have been already verified experimen- 
tally. However, many unexplained predictions remain 
at low temperatures both in the linear and the nonlin- 
ear regimes, as well as regarding the nonlinear regime 
at high temperatures. We hope that the rather precise 
derivations offered in this paper will guide such future ex- 
perimental work. Furthermore, some of the morphologies 
expected in the nonlinear regime need further theoretical 
understanding as well, allowing for the continuation of 
this inquiry. With the dramatic advances in computer 
speed, the understanding of some of these questions, ei- 
ther through numerical integration of the continuum the- 
ory or through discrete models, might be not too far. 
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APPENDIX A: 

The algebraic relation between the coordinates of the 
laboratory frame and the local frame, depicted in Fig. 
follows from the definitions given in point (i) of Sect. 
Accordingly, the unit vector along the Z axis is the 
normal at point O 



{-dxh,-dyh, 1) 
V9 



(Al) 



where g = 1 + {dxh)^ + (dyh)^. The vector m drawn on 
Fig. 1^ has components 

ifi = (sin 6*, 0, cos 9). 



(A2) 



Therefore, the unit vector along the Y axis reads: 

I X TO 1 



Y 



\n X TO 



{—{dyh) COS 0, sin ( 



'{dxh) cosd, (dyh) s'md), 



(A3) 



and finally (Al) and (A3) yield for the unit vector along 
the X axis: 

X ^Y X Z 
1 



(sine + (dxh) cos6'(9„/i)2sin6', 
gsimp 

(dyh) cos9 — (dxh)(dyh) sin 6*, 

(dxh) sme+((dxh)^ + (dyh)^) cose) . 



(A4) 



The matrix Ai defined in Eq. (|3^) and which relates the 
coordinates in the local and laboratory frames reads 
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where s — sin0, c — cose, and 
r=^(s + cdxh)^ + (dyh)^. 
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the Gaussian integrals in this formula we ob- 
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APPENDIX B: 



If we perform a small A„n 
obtain 



expansion in Eq. (pfl) we 



Taking into account Eqs. ( |33a| ) , (33b) relating the local 
((fi) and the global (e) angles of incidence through the sur- 
face slopes dxh, dyh, a small slope approximation leads 
to 
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which are the 6* ^ hmit of Eqs. (g3aD and ( p3b|) The 
small gradient expansion performed on Eq. (C2) now 
gives 



sin tp ~ ^ {dxhY + {dyhy. 



(C3) 



Using Eqs. (C3) in the expansions leading to Eq. ([40[), it 
can be seen that the expressions obtained for the coeffi- 
cients indeed are the 9 ^ Q limit of Eqs. (|4l|)-(|5F 
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r2o(e)-r("o^(e) + r«(0)(9./i), 
ro2(e)-r(°)(0) + r«(0)(9,.;i), 

+ aUyc'{s'-c')+alal 

X sV(7s2-5c2)+a;^c4(5s2-c2)} 

rg\9) = ^ {-2als' ~ ay + alal 

X (s2 + s2c2)_a4^4(^2^^2^2) 



I 42,'r-2 o2 2n 1^2421 

+ a^a^[os — Ss c ) + Aa^a^c | 



r(°)(e) = 



42^(^) 



2/ 
P ■ 



In the above expressions we used the notations 
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APPENDIX D: 



The solution corresponding to a rotated ripple struc- 
ture follows from Eq. (|83|). Indeed, in the absence of the 
S^x and ^j, terms, if we consider a solution of (^) of the 
form h{x,y,t) = g{x — vy,t) with v an arbitrary con- 
stant, the surface morphology evolution equation takes 
the form 



dtg 



-vq + -fdig +ivx +v Vy)di g 
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,.2\2q4. 



+ ^(A, + v^\y){digf + {Qi + v^n2)dfg 



-K{1 
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'd!g, 



(Dl) 



where g(l) — g{x — vy) is the steady wave solution [104|. 
From (Dl) it follows that the nonlinearity vanishes when 
Xx + v'^Xy = 0, or V = \J—\x/Xy. In this case we ob- 
tain an exponentially growing ripple structure with rip- 
ples forming an angle 9c = tan^^(w) = tan^^(^— A^/A.y) 
with respect to the x axis. 



n r 2 2 I 22 
c = cos 9, J = a^s + a^c . 
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APPENDIX C: 



Equations ( |34aD -( p4b| ) relating the incidence angle as 
measured in the local and laboratory reference frames ap- 
ply only in the off- normal incidence case (0 7^ 0). In the 
following we derive the correct expressions for the normal 
incidence case {9 = 0). Indeed, if = 0, the vectors n 
and TO shown in Fig. ^ are given by 



{-dxh,-dyh, 1) 
V9 



(0,0,1) 



(CI) 



Proceeding now as in (|33aD-(33b), we obtain 
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